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ABSTRACT 


This  thesis  presents  a  theoretical  investigation  of  the  coupled 
heat  transfer  problem  of  conduction  through,  and  laminar  free 
convection  from,  a  radial  triangular  fin  rotating  synchronously  with 
its  surroundings  at  a  constant  high  speed  about  a  horizontal  axis. 

When  the  Biot  number  is  assumed  to  be  small  and  the  Prandtl 
number  is  of  order  one,  two  essential  parameters  are  revealed  through 
an  order-of-magnitude  analysis.  These  are  the  ratio  of  the  fin 
length  to  the  distance  from  the  fin  tip  to  the  center  of  rotation  (y) ; 
and  the  product  of  the  fin  slenderness  ratio,  the  thermal  conductivity 
ratio  and  the  Rayleigh  number  to  the  fourth  power  (C) . 

The  overall  problem  was  solved  numerically  after  a  transformation 
favored  by  a  "local  similarity"  technique  had  been  applied.  The 
interfacial  temperature  and  heat  flux  between  solid  fin  and  the  fluid 
were  successfully  matched  to  a  satisfactory  degree. 

The  effects  of  parameters  C  and  Y  on  the  validity  of  similarity 
solutions,  Nusselt  number,  fin  effectiveness,  and  the  corresponding 
boundary  layer  profiles  were  studied  over  the  ranges  of  0.01  <  C  <  10 
and  0  <  y  <  0.3.  All  the  results  appear  to  be  interesting  and  reason¬ 


able  . 


I 


ACKNOWLEDGEMENTS 


The  author  wishes  to  express  the  sincere  appreciation  to 
the  following: 

Prof.  G.S.H.  Lock  for  his  patient  supervision  and 
assistance  during  the  course  of  this  thesis, 
the  National  Research  Council  for  financial  support 
under  grant  no.  A-1672. 

ray  wife,  Lily,  for  her  understanding  during  my  thesis 
writing  and  her  patience  in  typing  this  thesis. 


- 


V 


TABLE  OF  CONTENTS 


PAGE 

CHAPTER  I  INTRODUCTION  .  1 

CHAPTER  II  THEORETICAL  ANALYSIS  .  5 

2.1  BOUNDARY  LAYER  PROBLEM  .  5 

2.1- 1  MATHEMATICAL  FORMULATION  .  5 

2.1- 2  LOCAL  SIMILARITY  TRANSFORMATION  .  9 

2.2  CONDUCTION  PROBLEM  OF  FIN  . 10 

2.3  MATCHING  OF  THE  COUPLED  CONDUCTION  AND 

CONVECTION  SYSTEMS  .  14 

CHAPTER  III  RESULTS  AND  DISCUSSION  .  17 

3.1  SOLUTIONS  OF  CONVECTION  PROBLEM  .  17 

3.1- 1  DEPARTURE  FROM  SIMILARITY  . 17 

3.1- 2  EFFECT  OF  PARAMETER  y  . 19 

3.1- 3  EFFECT  OF  PARAMETER  C  . 21 

3.1- 4  NUSSELT  NUMBER  . 23 

3.1- 5  FIN  EFFECTIVENESS  e  . 27 

3.2  SOLUTIONS  OF  CONDUCTION  EQUATION  .  30 

CHAPTER  IV  CONCLUSIONS  AND  RECOMMENDATIONS  .  39 

4.1  CONCLUSIONS  . 39 

4 . 2  RECOMMENDATIONS 


42 


' 


vi 


PAGE 

REFERENCES  .  44 

APPENDIX  A  NORMALIZATION  .  46 

APPENDIX  B  LOCAL  SIMILARITY  TRANSFORMATION  .  52 

APPENDIX  C  QUASI-ONE-DIMENS IONAL  PROBLEM  OF  FIN  .  63 

APPENDIX  D  FORMULATION  OF  THE  LIMITING  CASE  y  =  0  . 69 

APPENDIX  E  NUMERICAL  TECHNIQUE  .  75 


' 


vii 


LIST  OF  ILLUSTRATIONS 

FIGURE  PAGE 

2.1  SYSTEM  OF  RADIAL  ROTATING  FIN  . 6 

2.2  DIMENSIONS  OF  TRIANGULAR  FIN  .  11 

3.1  DEPARTURE  FROM  SIMILARITY  . 18 

3.2  EFFECT  OF  y  .  20 

3.3  EFFECT  OF  C  .  22 

3.4  HEAT  TRANSFER  RELATION  WITH  y  AS  THE  PARAMETER .  25 

3.5  HEAT  TRANSFER  RELATION  WITH  C  AS  THE  PARAMETER  .  26 

3.6  FIN  EFFECTIVENESS  .  28 

3.7  EFFECT  OF  C  ON  INTERFACIAL  TEMPERATURE  .  31 

3.8  EFFECT  OF  y  ON  INTERFACIAL  TEMPERATURE  .  33 

3.9  RATE  OF  CONVERGENCE  .  35 

3.10  INTERRACIAL  TEMPERAUTRE  PROFILES  FOR  CHECKING  THE  POSSIBILITY 

OF  POWER  LAW  ASSUMPTION,  c}>0  =  xU  .  36 

3.11  INTERFACIAL  TEMPERATURE  PROFILES  FOR  CHECKING  THE  POSSIBILITY 

TTIX 

OF  EXPONENTIAL  ASSUMPTION,  <J>0  =  e  . 37 


, 


NOMENCLATURE 


£,x,x 

radial  distance  from  fin  tip 

p>pd 

pressure  departure 

R 

radial  distance  from  axis  of  rotation 

y»Y 

lateral  distance  from  fin  surface 

?,u 

local  similarity  variable 

u,U 

radial  velocity 

v,  V 

circumferential  velocity 

stream  function 

<M,9,T 

temperature 

O) 

angular  velocity 

A 

cross  sectional  area  of  fin 

L 

fin  length  (radially) 

W 

root  half-width  of  fin 

g 

gravitational  acceleration 

k 

thermal  conductivity 

$ 

thermal  expansion  coefficient 

V 

momentum  diffusivity 

K 

thermal  diffusivity 

C 

P 

constant  pressure  specific  heat 

h 

heat  transfer  coefficient 

D 

span  of  fin 

/ 

IX 


t ,  t  t ime 

a  angle  from  vertical  line  to  radial  axis  in  counterclockwise  direction 

Nu  Nusselt  number  (fiL/k^) 

2  3 

Ra  Rayleigh  number  ($to  R^G^L  /vk) 

Pr  Prandtl  number  (v/k) 

Bi  Biot  number  (hW/k  ) 

s 

2 

Ek  Ekman  number  (v/toL  ) 

2 

Fr  Froude  number  (R^co  /g) 

2 

Os  Ostrach  number  (3R  to  L/c  ) 

t  P 

y  length-radius  ratio  (L/R^) 

1/4 

C  dimensionless  group  (Ra  Lk  /Wk  ) 

JL  S 


SUBSCRIPTS 

0 

fin 

r,t 

fin  root,  tip 

00 

outside  boundary  layer 

co 

CH 

fluid,  solid 

y>* 

differentiation  with 

respect  to  y, 

i 

differentiation  with 

respect  to  n 

CHAPTER  I 


INTRODUCTION 


Free  convection  phenomena  arise  in  a  non-isothermal  fluid  because 
of  body  forces  which  generate  a  density  difference.  For  example,  when 
a  vertical  surface  generates  a  lateral  temperature  gradient  in  its 
surrounding  fluid,  a  boundary  layer  forms  adjacent  to  the  surface. 

Although  many  papers  have  dealt  with  steady  laminar  free  convection 
from  a  vertical  wall  subjected  to  a  parallel  uniform  body  force  field, 
only  those  solutions  with  similarity  characteristics  have  been  inten¬ 
sively  worked  out,  e.g.,  references  [1-5]*.  It  was  analyzed  in  Yang's 
work  [4]  that  similarity  solutions  are  possible  only  when  the  wall 
temperature  distribution  satisfies  some  prescribed  forms.  In  fact, 
similarity  solutions  are  very  special  cases  in  the  physical  world. 

When  the  wall  temperature  profile  does  not  exhibit  similarity  forms, 
the  partial  differential  equations  governing  the  boundary  layer  system 
can  no  longer  be  transformed  into  ordinary  differential  forms.  Hence 
non-similar  problems,  which  are  much  more  difficult  to  solve,  are 
encountered . 

Obviously,  instead  of  a  uniform  body  force  (such  as  gravity),  free 
convection  can  also  be  induced  by  centrifugal  and  Coriolis  forces. 
Although  Kreith  [6]  gave  a  detailed  review  of  many  papers  on  convective 
heat  transfer  in  rotating  systems,  very  few  publications  of  laminar  free 
convection  problems  on  a  vertical  surface  rotating  synchronously  with 
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surroundings  are  available.  In  1963  and  1964,  Lemlich  and  co-workers 
[7-9]  investigated  the  problem  of  an  isothermal  cold  plate  rotating 
with  its  leading  edge  at  the  center  of  rotation.  Meanwhile,  Catton  [10] 
applied  an  integral  technique  to  the  rotating  isothermal  plate  problem 
(with  large  Prandtl  number)  but  shifting  the  leading  edge  to  a  certain 
distance  from  the  center  of  rotation.  It  was  not  until  1969,  that  a 
similarity  solution  of  the  former  cold  plate  problem  was  presented  by 
Manoff  and  Lemlich  [11]. 

It  should  be  noticed  that  these  cold  rotating  plate  problems 
contain  a  tacit  assumption,  i.e.,  the  gravitational  field  does  not  exist. 
Without  this  assumption,  the  gravitational  force  will  periodically 
affect  the  system  in  the  vicinity  of  the  center  of  rotation,  a  horizontal 
axis,  and  no  leading  edge  of  a  steady  boundary  layer  can  take  place. 
Hence,  only  Catton’ s  attempt  of  shifting  the  leading  edge  from  the  center 
of  rotation  appears  justifiable  to  the  problems  under  the  simultaneous 
influence  of  gravitational  force. 

From  the  practical  viewpoint,  the  heat  transfer  from  a  hot  surface 
is  more  significant  than  from  a  cold  surface.  The  application  of 
thermal  fins  to  increase  the  cooling  rate  of  a  heated  body  is  a  good 
example.  In  a  uniform  body  force  field  parallel  to  the  vertical  surface, 
the  physical  phenomenon  of  a  hot  plate  can  easily  be  simulated  from  the 
study  of  a  cold  plate  (and  vise  versa)  by  choosing  the  opposite  end  of 
the  plate  as  the  leading  edge  and  reversing  the  longitudinal  axis. 
However,  when  the  body  force  is  non-uniform,  this  "reversibility"  can 
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not  hold  anymore.  As  an  example,  a  heated  radial  fin  spinning  together 
with  its  surroundings  about  a  horizontal  axis  in  a  gravitational  field 
is  somewhat  different  from  the  aforementioned  cold  rotating  plate 
problem.  Besides  the  fact  that  a  stable  leading  edge  can  now  form  at 
the  tip  of  the  thermal  fin,  two  substantial  differences  are: 

(1)  the  centrifugal  force  decreases  linearly  along  the  flow 
direction,  whereas  it  increases  in  the  cold  rotating  plate 
problem. 

(2)  the  similarity  assumption  of  the  cold  plate  temperature 
distribution  may  not  satisfy  the  surface  temperature  of  a 
rotating  fin  heated  at  the  root,  since  the  boundary  layer 
problem  is  always  coupled  with  the  conduction  problem  in  the 
fin. 

The  first  difference  between  these  two  problems  results  in  a 
radical  difference  in  the  equations  of  motion,  while  the  second  gives 
rise  to  a  problem  of  interfacial  temperature  and  heat  flux  matching. 

In  1968,  when  Lock  and  Gunn  [12]  presented  the  similarity  solution 
of  a  tapered  fin  parallel  with  uniform  body  force,  the  interfacial 
temperature  between  the  solid  and  liquid  was  matched  by  imposing  the 
similarity  assumption  of  the  boundary  layer  problem  upon  the  quasi-one- 
dimensional  conduction  equation  of  the  fin.  The  matched  solutions  were 
still  constrained  to  the  limited  case  of  power  law  temperature  distri¬ 
butions  along  the  fin.  A  coupled  heat  transfer  problem  without  power 
law  restriction  on  the  solution  was  studied  by  Kelleher  and  Yang  [13]. 


4 


In  that  work,  the  two-dimensional  conduction  problem  was  treated 
analytically  with  Fourier  transforms  and  the  boundary  layer  problem 
was  solved  with  Gortler-type  series.  Though  the  interfacial  tempera¬ 
ture  and  heat  flux  were  formally  matched  by  comparing  the  series 
solutions,  their  non-similarity  results  were  stated  to  be  valid  only 
near  the  leading  edge. 

In  the  present  analysis,  the  steady  coupled  heat  transfer 
problem  of  a  heated  radial  fin  rotating  simultaneously  with  the 
ambient  air  is  studied.  The  fin  is  located  at  a  certain  distance 
from  the  center  of  rotation  so  that  the  centrifugal  acceleration  can 
override  the  periodically  changing  effect  of  gravitational  force. 

One  can  verify  with  an  order  of  magnitude  that  the  Coriolis  force 
effect  along  the  flow  direction  becomes  negligible  when  the  rotating 
speed  is  high,  and  hence  the  centrifugal  force  dominates  the  overall 
region  of  interest. 

A  suitable  transformation  is  then  introduced  so  that  the  coupled 
governing  equations  can  be  solved  numerically  with  a  "local  similarity" 
technique  [17].  The  matching  problem  of  interfacial  temperature  and 
heat  flux  is  thoroughly  resolved  in  a  few  iterative  loops. 
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CHAPTER  II 


THEORETICAL  ANALYSIS 


2.1  BOUNDARY  LAYER  PROBLEM 


Laminar  free  convection  from  an  outward-projecting  fin  in  a 
centrifugal  force  field  is  studied  theoretically  for  an  arbitrary 
surface  temperature  increasing  monotonically  from  the  tip  to  the 
root . 

2.1-1  MATHEMATICAL  FORMULATION 


The  system  of  coordinates  is  displayed  in  figure  (2.1).  To  make 
the  problem  mathematically  tractable,  the  following  assumptions  were 
made: 

(1)  The  fluid  is  Newtonian. 

(2)  The  flow  is  laminar. 

(3)  The  fluid  is  mechanically  incompressible. 

(4)  The  flow  is  two  dimensional. 

(5)  Fluid  properties,  except  density  in  the  buoyancy  effect, 
are  constant. 

(6)  There  is  no  heat  generation. 

(7)  The  fluid  is  assumed  to  be  of  sufficient  extent  that  wall 
effects  are  insignificant. 

(8)  Radiation  is  neglected. 


CENTER  OF  ROTATION 
OF  SYSTEM 


Figure  2.1  System  of  Radial  Rotating  Fin 
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(9)  The  mechanical  work  against  the  body  force  field  is  assumed 
to  be  small. 

(10)  The  fin  slenderness  is  so  large  that  the  effect  of  inclination 
and  curvature  become  negligible. 


Through  the  processes  of  normalization  and  ordering,  it  has  been 
shown  in  appendix  A  that  for  air,  the  problem  becomes  quasi-steady  at 
high  rotating  speed.  By  referring  to  equations  (A. 16) - (A. 20)  in  appendix 
A,  the  governing  equations  in  dimensionless  form  are 


(2.1) 


(2.2) 


(2.3) 


(2.4) 


and 


3y  3y2 


where  X  =  -1  or  +1  for  the  trailing  and  leading  faces,  respectively. 
Equation  (2.3)  indicates  that  the  Coriolis  force  merely  induces  the 
lateral  pressure  gradient  in  this  problem.  Since  the  pressure  term  has 
been  suppressed  from  equation  (2.2),  one  can  find  the  pressure  distribu¬ 
tion  without  any  difficulty  once  the  profiles  of  u  and  <f>  are  at 
hand.  Hence  equation  (2.3)  can  be  ignored  in  this  problem  unless  one 
is  interested  in  the  pressure  profile.  Hereinafter,  the  problem  treated 
is  that  of  solving  the  other  three  partial  differential  equations  with 


the  associated  boundary  conditions: 
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x  =  0  ; 


<J>  =  0 


y  =  0  (x  >  0)  ; 


u  =  v  =  0 


<P  =  <Po(x) 


(2.5) 


y  =co  ; 


u  =  0 


(p  =  0 


Looking  at  the  overall  problem  of  a  fin  heated  at  its  root,  the 
surface  temperature  4>q(x)  is  normally  fixed.  There  is  no  point  in 
obtaining  similarity  solutions,  even  if  there  are  any,  of  the  boundary 
layer  problem  by  forcing  oneself  to  prescribe  a  specific  form  of  this 
interfacial  temperature.  In  other  words,  a  scheme  for  obtaining  more 
general  results  than  similarity  solutions  is  required. 

The  ’’local  similarity"  concept  [17]  suggests  that  the  whole 
region  can  be  discretized  along  the  longitudinal  direction  into  parallel 
levels.  At  each  level,  the  integration  is  carried  out  laterally  across 
the  boundary  layer  while  the  longitudinal  derivatives  in  the  governing 
equations  are  approximated  with  finite  differences.  In  this  way,  it  is 
not  necessary  to  specify  the  interfacial  temperature  <J>q(x)  in  any  re¬ 
stricted  form  of  x.  On  the  contrary,  cf>Q  can  be  an  arbitrary  set  of 
discrete  values  increasing  monotonically  from  the  tip.  The  flexibility 
of  these  discrete  values  is  very  advantageous  to  the  problem  of  matching 
when  the  iterative  method  is  applied  in  a  numerical  program  as  will  be 


seen  later. 
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2.1-2  LOCAL  SIMILARITY  TRANSFORMATION 


A  stream  function  \p(x,y )  is  defined  by  the  relations 


u  = 


3 1); 


v  =  - 


3iJj 


(2.6) 


3y  ,  3x  , 

and  equation  (2.1)  is  automatically  satisfied.  In  consequence  of 
removing  the  singularity  due  to  <j>o(0)  <j)(0,0),  see  appendix  B,  the 

"local  similarity"  concept  leads  to  the  following  transformation: 

5  -  yx  (2.7) 

T  1/4 


n  =  y 


d5 


3Y 


I 


f(5,n)  = 


3y 


G01/3(S)dS 


G01/3(S)dS 


-3/4 


<Kx,y) 


(2.8) 


(2.9) 


<K5,n)  =  G0  1(?)  4>(x,y) 


(2.10) 


where 


G0  (?)  =  <j>o  (x)  .  (2.11) 

As  outlined  in  appendix  B,  using  equation  (2.6)  and  the  "local 
similarity  transformation",  equations  (2.2)  and  (2.4)  yield 


Pr  (f,,f  +  $)  +  ff??  -  |  (f')2  =  -  ~ 


d5 


i 


dGo"1/3«>-  l  G01/3(S)dS 


? 

(f’) 


4 


•1/3 


a) 


1/3 


(S)dS 


f'  If"  "  f”  If)  +  Pr5$  ’  (2‘ 


and 

’  +  f$T  = 


-4 


dGn  1/3(P  f  „  1/3 


d5 


'0 


(S)  dS 


f 


Go  1/3(5)  r G01/3(S)dS 
~v 


85 


-  M) 

9  5/, 


12) 


(2.13) 


:«,•**  t®)  .  of|r,  ».f  W  (M  s**1 


where  the  primes  denote  partial  differentiations  with  respect  to  p. 

In  the  right-hand-sides  of  equations  (2.12)  and  (2.13),  each 

coefficient  is  a  function  of  £,  if  the  ^-derivatives  are  replaced  by 

finite  differences,  these  two  equations  become  ordinary  differential 

equations  for  each  specific  value  of  £  and  a  "local  similarity" 

problem  is  in  hand.  Provided  that  the  derivatives  of  f  and  $  are 
dG  "1/3 

finite  and  — - -  is  bounded  at  the  leading  edge,  the  right-hand-side 

of  equations  (2.12)  and  (2.13)  vanish  as  £  =  0.  With  the  associated 
boundary  conditions: 

p  =  0  ;  f=fT  =$-1=0  , 

n=oo  ;  f  ’  =  $  =  0  , 

there  is  no  difficulty  in  starting  the  integration.  The  solution  of 
these  equations  along  £  =  0  will  give  the  "initial"  values  of  this 
boundary  layer  problem. 

2 . 2  CONDUCTION  PROBLEM  OF  FIN 

Heat  transfer  within  a  fin  is  demonstrated  in  appendix  C  to  be 
a  quasi-one-dimensional  conduction  problem  when  the  Biot  number  (Bi) 
is  much  less  than  order  one. 

The  basic  notation  and  scheme  of  coordinates  of  this  triangular 
fin  problem  are  shown  in  figure  (2.2). 


§Jj 
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Figure  2.2  Dimensions  of  Triangular  Fin 
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In  addition  to  the  assumption  of  Bi  <<  0[1],  it  will  be  assumed 
that  there  is  no  heat  generation  in  the  fin  and  the  properties  of  the 
fin  are  constant.  From  the  order-of-magnitude  procedure  in  the  boundary 
layer  problem,  it  is  known  that  at  high  rotating  speeds,  the  temperature 
as  well  as  the  velocity  profiles  of  boundary  layer  flow  are  symmetric 
about  the  X-axis.  Therefore,  the  convective  heat  transfer  coefficients 
h^(X)  and  h^CX),  shown  in  figure  (2.2),  are  equal  and  can  be  represented 


by 


_K  9T(X,0) 
h! (X)  =  h2(X)  =  T0(X)  -  T 


(2.14) 


The  temperature  Tq(X)  of  the  triangular  fin  satisfies  the  equation 


d2Tp (X)  dTp (X)  ,  L^f  9T(X,0)  = 

,v 2  dX  W  K  3Y 

dX  s 


=  0 


(2.15) 


Normalizing  this  equation  with  those  characteristic  values  obtained  in 
appendix  A  and  introducing  the  "local  similarity  transformation"  yield 
a  nonlinear  ordinary  differential  equation,  i.e. 


d2Gn(Q 

d£2 


,  dGn(g) 

d? 


f  \  3/4 

C$’(?,0)  G01/3(S)dsj 


G0(5)  =  0 


(2.16) 


where  $T(C,0)  is  one  of  the  eigenvalues  from  the  solution  of  the  boun¬ 
dary  value  problem  and  C  is  a  parameter,  defined  [12]  as 

c  -  k  JE  Ra1/A 

C  -  W  K  Ra 
s 


(2.17) 


’ 
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It  is  clear  that  one  of  the  boundary  conditions  for  this  second- 


order  differential  equation  (2.16)  is 

g0(y)  =  1 

at  the  root  of  the  fin.  The  other  boundary  condition  needs  more  con¬ 
sideration.  It  is  shown  in  appendix  C  that  when  the  slenderness  ratio 

is  much  larger  than  order  one,  either 

w 

G0(O)  =0  or  dG^(0)-  =  0 

is  a  good  approximation  to  the  boundary  condition  at  the  tip.  In  view 
of  the  generality  of  the  problem,  it  is  undesirable  to  be  constrained 
to  the  case  of  a  very  thin  fin.  A  boundary  condition  which  is  also 
realistic  to  a  triangular  fin  with  smaller  slenderness  is  required. 
Equation  (2.16)  can  be  rewritten  as 


(2. 


and  it  can  be  regarded  as  a  first-order  nonlinear  equation  in  £ 


dGa 

d£  ’ 


and  therefore  one  boundary  condition  on  this  variable  is  sufficient. 
dC 

Since  g  indicates  the  total  amount  of  heat  conducted  across  the 
at, 

fin  section  at  a  distance  E,  from  the  tip,  it  tends  to  zero  as  the  area 

dG 

of  the  cross-section  becomes  zero  at  the  tip.  Hence  =  0  at  the 

tip  can  be  the  boundary  condition  of  equation  (2.18).  Once  the  solu- 
tion  of  is  obtained,  G0  =  1  can  be  used  as  the  boundary  condition 

dq 

at  the  root  for  further  integration  in  obtaining  Gq. 


t 
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2.3  MATCHING  OF  THE  COUPLED  CONDUCTION  AND  CONVECTION  SYSTEMS 

The  whole  problem  has  been  shown  to  be  governed  by  two  boundary 
layer  equations,  (2.12)  and  (2.13),  and  one  conduction  equation,  (2.18). 
All  these  three  governing  equations  are  mutually  coupled.  The  inter¬ 
facial  heat  flux  between  solid  fin  and  fluid  has  been  expressed  in 
equation  (2.14)  and  in  succession  equations  (2.16)  and  (2.18)  were 
formed.  Hence  the  only  interfacial  matching  left  to  be  done  is  on 
the  temperature,  i.e.,  the  assumed  wall  temperature  Gq(5)  of  the  boun¬ 
dary  layer  problem  must  satisfy  the  conduction  equation  of  the  fin 
and  its  associated  boundary  conditions. 

As  was  stated  before,  Go (O  is  a  fixed  but  unknown  wall  tempera¬ 
ture,  and  therefore  a  logical  way  of  solving  this  problem  is  to  dis¬ 
cretize  the  ^-coordinate  into  certain  levels  and  assume  a  set  of  values 
for  Gq(5)  at  these  levels.  These  assumed  values  decrease  monotonically 
from  Gq (y)  =  1  to  Gq(0)  ^  0.  By  means  of  this  initial  guess,  the  three 
aforementioned  coupled  equations  can  be  solved  to  give  a  new  set  of 
values  for  Gq  (?) •  After  comparison  and  modification,  convergence  of 
this  interfacial  temperature  G0(?)  to  a  unique  distribution  can  be 
procured  through  an  iterative  method. 

The  numerical  technique  is  too  lengthy  to  be  presented  in  this 
section  and  its  details  are  described  in  appendix  E:  only  the  outlined 
procedures  are  given  in  sequence: 


, 


(1)  With  the  initial  set  of  Gq(5)  as  an  input  to  the  boundary 
layer  problem,  equations  (2.12)  and  (2.13)  can  be  solved 
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(2) 


simultaneously  at  each  £-level  where  finite  difference 


approximations  replace  the  ^-derivatives .  The  integrations 
in  the  q-direction  are  carried  out  with  a  fourth-order 
Runge-Kutta  method  and  convergence  is  obtained  by  means 
of  the  technique  developed  by  Nachtsheim  and  Swigert  [14]. 
Once  the  boundary  layer  problem  has  been  solved,  the 
eiganvalues  $’(£,0)  at  each  level  are  known  and  can  be 
substituted  into  equation  (2.18).  The  bracketed  group  in 
this  equation  represents  the  convective  heat  transfer 
coefficient  h(GQ,£).  In  order  to  circumvent  the  non¬ 
linearity  of  this  equation,  the  assumed  Gq(£)  distribution 
is  used  in  calculating  this  coefficient.  A  fourth-order 


Runge-Kutta  method  is  then  employed  to  integrate  equation 

(2.18)  from  the  root  to  give  Gq(£)  and  g;-  -0-  at  each  level 

except  the  leading  edge,  because  of  a  singularity.  This 
latter  difficulty  is  overcome  when  Gq(£)  is  obtained  by 
integrating  equation  (2.16)  simultaneously  with  equation 

(2.18) .  The  integration  starts  with  the  boundary  condition 

dG, 


Gq  (y)  =  1  and  a  guessed  value  of  £- 


which  is  gradually 


5=y 


rQvy,  -  —  ^  ~~  Sd£ 

modified  by  means  of  a  shooting  method  until  the  other 

dG, 


boundary  condition  £ 


error  less  than  10 


-5 


=  0  is  satisfied  within  an 


5=0 


- 
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(3)  A  new  Gq(5)  is  now  obtained  by  halving  the  two  most  recent 
solutions  and  is  substituted  into  the  bracketed  groups  of 
equations  (2.16)  and  (2.18)  to  calculate  a  new  set  of 
h(Go,0.  The  iteration  continues  until  Gq(£)  has  converged 
satisfactorily  for  the  same  input  values  of  $'(£>0). 

(4)  This  converged  Gq (0  is  regarded  as  a  better  input  to  the 
boundary  layer  problem  and  the  whole  procedures  starts  all 
over  again.  Iteration  between  the  conduction  and  convection 
systems  stops  when  the  normalized  interfacial  temperature 

G0 (O  has  converged  at  every  point  to  within  an  error  of 
less  than  10  Since  the  temperature  distribution  Gq(£) 

is  obtained  and  the  interfacial  heat  flux  is  matched  auto¬ 
matically,  the  problem  is  completely  solved. 


I 


CHAPTER  III 


RESULTS  AND  DISCUSSION 

Throughout  the  study,  all  the  calculations  are  based  on  taking 
air  as  the  fluid  and  0.72  as  its  Prandtl  number.  By  referring  to  the 
governing  equations  (2.12),  (2.13),  (2.16)  or  (2.18)  and  their  associ¬ 
ated  boundary  conditions,  it  is  evident  that  C  and  y  are  the  principal 
parameters  of  this  rotating  radial  fin  problem.  In  this  study,  solu¬ 
tions  are  developed  with  four  values  of  y  (0,  0.1,  0.2  and  0.3)  and 
six  values  of  C  (0.01,  0.1,  1.0,  3.0,  5.0  and  10.0).  Matched  results 
can  be  separated  into  two  parts:  convection  and  conduction. 

3.1  SOLUTIONS  OF  CONVECTION  PROBLEM 

In  a  boundary  layer  problem,  interest  most  often  rests  in  the 
velocity  and  temperature  profiles  of  the  fluid.  In  addition,  the 
Nusselt  number  and  fin  effectiveness  are  also  worth  investigating 
from  a  heat  transfer  point-of-view. 

3.1-1  DEPARTURE  FROM  SIMILARITY 

Temperature  and  velocity  profiles  at  three  E,  levels  along  a 
rotating  fin  with  C  =  1.0  and  y  =  0.3  are  shown  in  figure  (3.1).  If 
similarity  existed,  all  the  velocity  profiles  would  overlap,  as  would 
the  temperature  profiles.  However,  figure  (3.1)  shows  that  the  depar¬ 
ture  from  similarity  is  quite  significant  in  velocity  but  not  in  tem¬ 


perature  . 


Hence,  considering  the  heat  transfer  rate  from  a  rotating  fin 
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Figure  3.1  Departure  from  Similarity 
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with  small  y,  the  similarity  solution  obtained  locally  at  the  leading 
edge  provides  a  good  approximation.  This  may  not  be  suitable  when  y 
is  large,  since  the  variation  between  temperature  profiles  appears  to 
increase  as  E,  (through  y)  increases. 

3.1-2  EFFECT  OF  PARAMETER  y 


Figure  (3.2)  shows  the  velocity  and  temperature  profiles  at  the 
root  of  the  fin,  i.e.,  x  =  1.0,  for  four  values  of  y  with  C  =  1.0. 

It  is  obvious  that  the  temperature  profiles  are  less  sensitive  than 
velocity  profiles  to  the  variation  of  y.  Attention  must  be  paid  to 
the  limiting  case  of  y  =  0.  As  y  tends  to  zero,  the  body  force 
approaches  a  constant  value  and  the  problem  becomes  that  of  a  static 
fin.  In  other  words,  the  case  of  y  =  0  simulates  a  downward-projecting 
fin  in  a  gravitational  force  field  [12].  Appendix  D  presents  the  for¬ 
mulation  of  the  governing  equations  when  y  =  0,  and  shows  that  non¬ 
similarity  normally  exists  for  this  static  fin  problem.  Only  in  the 
very  special  case  when  C  =  0  and  an  isothermal  surface  is  formed  is 
the  similarity  solution  possible.  Since  C  =  1.0  is  chosen  in  this 
figure,  the  presented  profiles  with  y  =  0  are  the  local  similarity 
solutions . 

Again,  the  temperature  profiles  may  deviate  significantly  from 
those  with  y  =  0  when  the  value  of  y  approaches  unity;  likewise  for  the 
velocity  profiles. 
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C  =  1.0 


C  =  Y  (x  =  1.0) 
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Figure  3.2  Effect  of  y 
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3.1-3  EFFECT  OF  PARAMETER  C 

Tills  effect  has  been  investigated  over  three  orders-of-magnitude 
for  C.  Salient  changes  in  both  velocity  and  temperature  profiles  at 
the  root  of  the  same  fin  are  shown  in  figure  (3.3),  where  y=0.3. 

It  can  be  seen  that  as  the  magnitude  of  C  gets  smaller,  both  profiles 
approach  asymptotically  to  those  for  C=0:  this  evidently  occurs  in 
the  vicinity  of  the  presented  profiles  with  C-0.01.  Such  asymptotic 
behavior  is  very  useful  in  practice  because  whenever  C  is  small,  the 
result  for  C^O  is  a  good  approximation.  It  is  shown  at  the  end  of 
appendix  D  that  the  fin  temperature  Gq(0=1  when  C=0,  and  the  matching 
problem  no  longer  exists.  Certainly,  such  a  limiting  solution  for 
free  convection  from  an  isothermal  surface  takes  less  effort  and  time 
to  attain,  and  is  a  convenient  approximation  when  C  is  small. 

It  is  worth  noting  again  that  the  magnitude  of  parameter  C 
depends  on  the  fin  slenderness,  the  conductivity  ratio  of  fluid  to 
solid  and  the  Rayleigh  number  defined  by 

3  R  w2  0  L3 

t  r 

Ra  =  - 

V  K 

Clearly,  the  variation  of  C  may  represent  many  possible  situations. 
However,  if  all  but  one  of  the  component  parameters  of  C  are  regarded 
as  fixed,  the  effect  of  each  component  parameter  can  be  investigated 
separately  through  the  effect  of  C.  For  example,  it  appears  in  figure 
(3.3)  that  the  temperature  profiles  in  the  vicinity  of  the  wall  are 
more  declivous  for  larger  values  of  C.  This  behavior  implies  that  the 


, 
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Figure  3.3  Effect  of  C 
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heat  transfer  can  be  increased  by  means  of  a  longer  fin,  or  with  a 
good  conductive  rluid.  In  addition,  this  heat  transfer  can  also 
increase  with  the  temperature  difference  between  the  root  and  the 
surrounding  fluid,  or  with  the  rotating  speed.  As  C  increases,  it 
appears  that  the  temperature  boundary  layer  becomes  thinner  while 
the  momentum  boundary  layer  becomes  thicker  and  the  velocity  profile 
flattens  out. 

3 . 1-4  NUSSELT  NUMBER 

It  is  well  known  that  when  free  convection  alone  is  considered, 
the  Nusselt  number  for  a  vertical  plate  problem  in  a  uniform  gravi¬ 
tational  field  is  a  function  of  Grashof  number  and  Prandtl  number, 
e.g.,  reference  [1],  In  this  study,  consideration  is  extended  to  the 
interfacial  temperature  and  heat  flux  matching  problem  in  a  centrifugal 
force  field.  It  is  reasonable  to  expect  that  the  Nusselt  number  will 
be  a  function  of  additional  parameters,  which  enter  into  the  problem 
as  a  result  of  either  the  inclusion  of  conduction  or  the  introduction 
of  rotation.  From  the  analysis,  it  is  clear  that  C  and  y  are  the 
additional  parameters.  Alternatively,  the  local  Nusselt  number  at  the 
root  of  the  fin  is  defined  as 


=  h  (L)  L 

h 


-Ra 


1/4 


9$(y  >  0) 


(3.1) 


Since  in  this  study, 


are  functions  of  y  and  C 


■■■  *  asH 
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for  air,  one  can  assume  that  Nu/Ra  ^  is  a  function  of  y  and  C  as 
well,  i.e., 

Nu/Ra  =  F(C, y) . 

The  parameter  Nu/Ra^^  is  given  as  a  function  of  C  for  various 
values  of  y  in  figure  (3.4).  Just  as  predicted  in  section  3.1-4,  the 
increase  of  C  due  to  the  variations  of  its  "sub-parameters"  can 
improve  the  heat  transfer  and  hence  increase  the  Nusselt  number. 

A  more  interesting  heat  transfer  relation  appears  in  figure  (3.5) 
where  Nu/Ra"^^  is  plotted  against  y  at  six  different  values  of  C. 
Least-squares  lines,  with  a  maximum  error  less  than  1%,  suggest  that 
Nu/Ra^^  changes  linearly  with  the  parameter  y.  Therefore,  Nu/Ra"*"^ 
can  be  expressed  by 

— =  F1(C)-F2(C)y,  (3.2) 

Ra 

where  functions  F^  and  F2  increase  with  the  parameter  C.  The  magnitudes 
of  F^  and  F2  at  six  C  levels  are  presented  in  table  below: 


c 

F1 

F2 

0.01 

0.389 

0.116 

0.10 

0.401 

0.126 

1.00 

0.484 

0.178 

3.00 

0.573 

0.205 

5.00 

0.623 

0.214 

10.00 

0.696 

0.228 

TABLE  I 
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Figure  3.4  Heat  Transfer  Relation  with  y  as  the  Parameter 
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Figure  3.5  Heat  Transfer  Relation  with  C  as  the  Parameter 
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Though  values  of  Nu/Ra  ^  for  y>0.3  were  not  obtained,  one  can 
confidently  extrapolate  the  curves  in  figure  (3.5):  this  would  not 
have  been  said  of  figure  (3.4).  How  far  they  can  be  extrapolated  is 
certainly  unknown,  but  they  can  not  go  beyond  y=l.  Considering  this 
limiting  case  of  y  approaching  unity,  the  extrapolations  using  table  I 
and  equation  (3.2)  always  indicate  non-negative  heat  transfer  which  is 
physically  consistent. 


3.1-5  FIN  EFFECTIVENESS  e 


The  fin  effectiveness  for  a  heat  transfer  coefficient  h  varying 


along  the  fin  surface  has  been  defined  in  reference  [12]  as 

L 


£  = 


h (X) 0O (X) dX 


0 


h  (L)  0  L 
r 


(3.3) 


To  express  e  in  terms  of  the  normalized  and  transformed  functions  and 
variables,  use  is  made  of  equations  (2.7),  (2.8),  (2.10)  and  (2.11). 
Thus 


r  J  i 

1/4  J 

r  ,-5  i 

^Jg01/3(5)  d5 

/  g02(O$'(?,0)  ~ 

3 -[  G01/3(S)dS 

Jo 

n  3/4 


d? 


$’(Y,0) 


(3.4) 


Calculated  values  of  £  for  different  ratios  y  are  plotted  in 
figure  (3.6)  as  a  function  of  C.  The  approximate  solution  of 
reference  [12]  is  also  presented  as  a  comparison  to  the  static  fin 
case  when  y->0.  It  should  be  noted  that  no  curve  has  intersected  the 
coordinate  of  C=0 :  the  minimum  value  of  C  used  for  each  solid  line 
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Figure  3.6  Fin  Effectiveness 
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was  0.01.  However,  the  fin  effectiveness  for  C=0  can  be  predicted 
if  the  curves  are  extrapolated  to  the  ordinate.  The  value  of  e  with 


C=0  is  the  fin  effectiveness  of  an  isothermal  fin  with  Kf=  oo .  This 

value  of  e  does  not  apply  to  the  case  where  K  =0  or  0  =0,  since 

s  r 


both  the  numerator  and  denominator  of  definition  (3.3)  vanish 
(appendix  D) ;  hence  e  is  undefined. 

In  figure  (3.6),  the  effect  of  G  is  again  shown  to  be  greater 
than  that  of  y.  At  first  glance,  it  may  seem  surprising  that  the 
fin  effectiveness  can  be  greater  than  unity  when  C  is  small.  After 
the  definition  of  e  is  investigated  more  closely,  it  is  not  unreason¬ 
able  to  expect  such  a  result.  When  C  is  small,  the  fin  surface 
approaches  an  isothermal  condition.  In  the  definition  (3.3),  h(L)  is 
no  longer  the  maximum  value  of  h(X)  and  the  ratio  £  can  certainly  be 
greater  than  unity.  More  convincing  arguments  can  be  obtained  from 
the  well  known  isothermal  surface  problem  in  a  uniform  gravitational 
field  which  simulates  the  case  of  Y=0  when  C=0 .  By  means  of  the 
ordinary  similarity  transformation  for  an  isothermal  surface  [1, 3,4,5], 
£  can  be  calculated  from  definition  (3.3): 


1 


Apparently,  this  is  another  check  on  the  present  study,  since 
when  the  curve  of  y=0  extends  to  the  C=0  axis,  it  just  meets  this 
analytically  obtained  value  of  e.  As  a  matter  of  fact,  the  fin 


effectiveness  £  of  reference  [12]  is 


’ 

» 

* 


30 


4 

£  5n-f3 

where  n  is  the  index  of  the  power  law  surface  temperature  distribution. 
An  isothermal  fin  (C  =  n  =  0)  thus  yields 

£  =  4/3  . 

That  is,  in  this  isothermal  case,  the  predictions  of  both  studies 
become  coincident. 

3.2  SOLUTIONS  OF  CONDUCTION  EQUATION 

As  stated  before,  the  interfacial  heat  flux  between  solid  and 
fluid  has  been  matched  automatically  in  the  formulation  of  the  fin 
equation.  Hence,  only  the  interfacial  temperature  needs  attention 
and  investigation.  Since  the  fin  equation  is  of  quasi-one-dimensional 
form  as  Bi<<0  [1],  its  solutions  are,  in  fact,  the  matched  interfacial 
temperature  profiles. 

Figure  (3.7)  shows  the  salient  effect  of  the  parameter  C  on  the 
matched  surface  temperature  profile.  One  important  point  must  be 
mentioned;  as  C  approaches  zero,  the  normalized  wall  temperature  <j)Q 
becomes  close  to  unity,  this  phenomenon  strongly  suggests  that  <J>o  =  1 
when  C  =  0.  That  is,  the  numerical  results  are  once  more  supported  by 
analysis.  Another  feature  worth  pointing  out  is  that  the  surface 
temperature  and  its  slope  at  the  leading  edge  approach  zero  when  C  is 
large.  Therefore,  the  ordinary  boundary  condition  at  the  fin  tip 

M°)  =  0  °r  “dx~  =  0 
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Figure  3.7  Effect  of  C  on  Interfacial  Temperature 
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is  a  good  approximation  in  the  case  of  large  C  value:  this  is  consistent 
with  the  predictions  of  appendix  C,  i.e.,  as  —  «0[1],  <j>0(0)=0  and 

JLi 

^  =  0.  The  consistency  can  be  proved  as  follows.  The  validity 

of  the  quasi-one-dimens ional  fin  equation  is  assumed  when 


_  .  hW  _  r ,  -i 
Bi  =  —  «  0[1]. 

s 

Introduction  of  the  normalized  variables  of  appendix  A  gives 

W  ^ 

(-)  C  «  0[1]. 

Now,  the  result  shows  that  tf>o(0)=0  and  =  0  can  be  satisfied 

whenever  C  >>  0[1],  This  implies  that  ^  <<  0[1].  In  fact,  j-  <<  0[1] 
is  the  necessary  condition  of  the  homogeneous  boundary  conditions 
<J>0(0)=0  and  =  0  for  a  conduction  problem  while  C  >>  0[1]  is 

for  a  coupled  problem  of  conduction  and  laminar  free  convection. 

Anyhow,  in  this  problem  if  one  intends  to  use  the  aforementioned 
homogeneous  boundary  condition,  correct  results  can  only  be  guaranteed 
when  C  is  at  least  two  orders-of-magnitude  greater  than  unity. 

In  order  to  study  the  effect  of  y  on  these  interfacial  tempera¬ 
tures,  the  profiles,  based  on  the  normalized  variables  x,  for  four 
different  y  values  are  presented  in  figure  (3.8)  at  four  C  levels. 

The  profiles  at  each  C  level  virtually  coinside,  and  the  maximum 
difference  among  them  is  less  than  1%.  This  is  a  very  interesting 
feature  from  the  economic  viewpoint,  since  even  if  the  difference  is 
not  from  the  calculating  errors,  any  normalized  interfacial  temperature 
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Figure  3.8  Effect  of  y  on  Interfacial  Temperature.  Each 
Overlapped  Curve  Includes  Four  Curves  for 
y  =  0,  0.1,  0.2,  0.3  . 
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for  small  y  can  be  considered  as  a  good  initial  guessed  profile  for 
higher  values  of  y,  for  the  same  value  of  C.  Experience  shows  that 
the  larger  the  value  of  y,  the  longer  the  computing  time  required; 
and  the  closer  the  initial  guessed  profile  to  the  finally  matched 
profile,  the  faster  the  convergent  rate.  Hence  in  such  cases,  a 
good  initial  guess  can  definitely  cut  down  the  computation  required 
for  a  large  y  value  quite  significantly. 

The  convergence  rate  of  the  interfacial  temperature  is  indicated 
in  figure  (3.9)  where  two  significantly  different  initial  guesses.  Ml 
and  Nl,  were  tested  at  C  =  1.0  and  y  =  0.2.  After  one  loop  of  matching 
between  the  conduction  and  convection  system,  they  gave  very  good 
approximations  (M2  and  N2  respectively)  to  the  final  converged  profile. 
This  figure  not  only  reveals  the  fast  matching  ability  of  the  technique 
used  in  this  study,  but  also  gives  more  confidence  in  convergence  from 
"wild"  initial  guesses. 

Finally,  the  possibility  of  approximating  the  similarity  solution 

to  the  vertical  surface  convection  problem  is  investigated  at  y  =  0.3 

in  figures  (3.10)  and  (3.11).  Both  graphs  suggest  that  the  interfacial 

temperature  profile  4>q(x)  matched  between  conduction  and  convection 

generally  bear  neither  the  form  of  xn  nor  e™X.  An  exception  to  this 

occurs  when  C  is  small,  where  cf> 0  is  close  to  unity  along  the  fin.  Again, 

from  a  practical  standpoint,  the  solution  for  an  isothermal  surface 

is  a  good  approximation  to  that  for  small  value  of  C.  At  least  it  can 
/ 

be  a  good  initial  guess  to  the  problem  when  C  is  small. 
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Figure  3.9  Rate  of  Convergence 
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Figure  3.10  Interfacial  Temperature  Profiles  for  Checking 

the  Possibility  of  Power  Law  Assumption,  <J>q  =  x11 
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C  =  0.01 


Figure  3.11  Interfacial  Temperature  Profiles  for  Checking  the 


mx 

e 


Possibility  of  Exponential  Assumption,  <J>q  - 
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Figure  (3.7)  reveals  that  <j>o(0)  tends  to  zero  as  the  value  of 
C  increases.  This  suggests  that  a  power  law  distribution  of  surface 
temperature  might  be  a  good  approximation  at  large  values  of  C. 

This  is  borne  out  by  figure  (3.11)  which  reveals  that  when  the  value 
of  C  is  much  greater  than  unity,  the  exponential  surface  temperature 
distribution  emX,  which  is  the  asymptote  for  xn  when  n  ->  00  [18],  is 
a  good  approximation  near  the  root. 


CHAPTER  IV 


CONCLUSIONS  AND  RECOMMENDATIONS 


4.1  CONCLUSIONS 


The  coupled  problem  of  quasi-one-dimensional  conduction  and 

non-similarity  convection  for  a  heated  triangular  fin  rotating 

together  with  surrounding  air  has  been  solved  successfully  and 

efficiently.  The  calculated  functions  and  their  derivatives  being 

of  order  one  indicates  that  the  normalization  procedure  is  valid. 

Although  the  simplified  boundary  layer  equations  (2.1)~(2.4),  obtained 

by  means  of  the  order  of  magnitude  method,  require  that  all  the  dimen- 

-1/4  -1 

s ionless  quantities  Ra  ,  Fr  ,  SQr>  and  Os  should  be  much  less  than 
unity,  a  realistic  example  in  appendix  A  has  demonstrated  that  there 
is  no  conflict  among  these  requirements. 

The  normalized  governing  equations  indicate  that  when  the 
angular  velocity  of  the  rotating  system  is  high,  the  periodic  effect 
of  the  gravitational  force  can  be  suppressed  and  the  system  can  reach 
a  quasi-static  state.  The  Coriolis  force  is  insignificant  in  the 
longitudinal  direction  along  the  fin,  and  only  contributes  to  the 
pressure  gradient  in  the  lateral  direction.  This  favorable  feature 
not  only  simplifies  the  mathematical  problem  but  also  indicates 
symmetry  of  the  boundary  layer  system  about  the  fin. 


I 
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A  local  similarity  transformation  with,  its  coefficients  being 
either  constant  or  functions  of  the  unprescribed  interfacial  tempera¬ 
ture  distribution  has  been  developed  so  that  the  interface  is  free 
from  the  ordinary  similarity  restriction  on  the  temperature  profile; 
e.g.,  a  power  law  assumption.  This  flexibility  of  the  interfacial 
temperature  makes  the  whole  transformation  suitable  to  the  problem  of 
interfacial  matching.  Analysis  shows  that  only  in  the  special  case 
of  an  isothermal  fin  in  a  uniform  gravitational  field,  can  the  tradi¬ 
tional  similarity  transformation  be  applied  to  such  a  problem.  Both 
transformations  exhibiting  the  same  forms  in  their  coefficients  can  be 
regarded  as  a  good  check  on  this  local  similarity  transformation. 

The  symmetry  of  the  boundary  layer  systems  on  both  sides  of  the 
rotating  fin  with  large  angular  velocity  permits  the  use  of  the 
quasi-one-dimensional  equation  for  the  conducting  fin.  It  has  been 
demonstrated  in  this  study  that  the  condition  for  a  quasi-one-dimensional 
fin  equation  is  Bi  <<  0[1].  Again,  an  example  shows  that  there  is  no 
conflict  between  this  condition  and  any  of  the  aforementioned  four  con¬ 
ditions  in  the  convection  problem.  A  new  type  of  boundary  condition 

with  more  generality  for  a  tapered  fin  has  been  introduced  at  the  tip 

dT  (0) 

of  the  fin.  The  traditional  boundary  condition,  T  (0)  =  0  or  — -r -  =  0, 

w  ux 

is  valid  only  when  the  parameter  C  of  this  problem  is  much  greater  than 
order  one. 

Truncation  errors  have  been  carefully  controlled  in  the  numerical 
approximations.  After  a  preliminary  test,  the  step  size  for  integration 
in  both  the  lateral  and  longitudinal  directions  were  chosen,  respectively, 
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to  be  0.04  and  0.02,  such  that  the  eiganvalues,  $'(5,0)  and  f”(?,0), 
of  the  boundary  layer  problem  were  within  an  error  of  1%  of  their  asymp¬ 
totic  values,  and  with  the  computing  time  within  a  tolerable  scale. 

In  the  beginning,  no  guarantee  could  be  made  on  the  possibility 
and  rate  of  convergence  of  interfacial  temperature.  More  confidence 
was  gained  after  several  trials,  as  in  figure  (3.9),  were  conducted. 

When  double  precision  is  used,  temperature  can  converge  to  consistency 
in  five  decimal  places  within  30  minute  computation  on  an  IBM  360/67 
machine.  The  interfacial  temperatures  presented  in  this  work  have 
converged  to  sixteen  decimal  places.  As  a  check  on  the  numerical 
program,  the  total  heat  flux  across  the  fin  surfaces  was  compared  with 
the  heat  flux  applied  at  the  root  of  the  fin  after  each  converged 
result  was  obtained.  Comparisons  of  these  values  showed  that  the 
maximum  discrepancy  is  of  the  order  of  1%. 

In  addition  to  Prandtl  number,  two  parameters,  y  and  C,  appear 
in  this  study.  The  effect  of  the  parameter  y  was  found  to  be  smaller 
than  that  of  C.  When  the  parameter  y  approaches  zero,  the  problem 
becomes  a  static  downward-proj ecting  fin  in  a  uniform  gravitational 
field.  Analysis  has  shown  that  the  similarity  solution  only  exists  in 
the  isothermal  static  fin  problem,  that  is,  when  C  =  0 . 

In  fact,  the  isothermal  fin  result  with  C  =  0  can  also  be  extended 
to  the  problem  in  a  centrifugal  force  field  (y  >  0) .  The  solutions  of 
this  work  indicate  a  favorable  feature  for  practical  purposes,  i.e., 
when  C  is  small,  solutions  with  C  =  0  can  represent  a  good  approximation. 
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As  one  can  see,  the  problem  of  matching  no  longer  exists  in  the 
isothermal  case  C  =  0  and  the  solution  is  much  easier  to  obtain. 

On  studying  the  effects  of  y  and  C  on  the  local  Nusselt  number 
at  the  root  of  the  fin,  a  linear  relation  between  Nu  and  y  was  obtained. 

717?  ■  Vc)  -  Vc>  *  • 

Ra 

The  fin  effectiveness  as  a  function  of  parameter  C  were  calculated 
at  four  y  values.  The  results  for  y  =  0  were  compared  with  the  solutions 
of  reference  [12].  If  the  calculated  fin  effectiveness  for  a  static  fin 
(y  =  0)  is  extended  to  C  =  0  from  the  range  of  C  in  this  study,  it  will 
match  the  analytical  result  for  an  isothermal  static  fin,  i.e.,  e  =  4/3. 
This  reveals  that  the  solution  of  this  work  meets  the  solution  of  re¬ 
ference  [12]  in  the  isothermal  case. 

4 . 2  RECOMMENDATIONS 

Firstly,  experimental  work  is  strongly  suggested  so  that  the 
results  of  this  work  can  be  compared  with  the  physical  evidence. 

If  time  of  computation  were  not  a  factor,  it  is  worth  trying  the 
problem  with  y  larger  than  0.3.  The  phenomenon  in  figure  (3.8),  that 
the  normalized  interfacial  temperatures  are  almost  independent  of  y 
over  its  range  of  investigation,  suggests  that  it  can  be  a  good  initial 
guess  to  the  problem  for  large  y,  and  the  computing  time  on  the  matching 
can  be  reduced  accordingly. 

The  order-of-magnitude  procedure  in  this  study  refers  to  a  fluid 
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with  Prandtl  number  of  order  one  and  in  the  calculations,  air  with 
Pr  =  0.72  has  been  considered  as  the  fluid.  One  can  carry  out  the 
same  calculation  for  different  fluids.  If  the  Prandtl  number  is  not 
of  order  one,  a  similar  order-of-magnitude  may  be  required,  and 
different  governing  equations  for  the  boundary  layer  problem  may  be 
formulated . 

Further  studies  can  be  made  either  on  the  ultra-high  rotating 
speed  problem  or  in  the  problem  with  small  angular  velocity.  The 
former  will  bring  the  heat  dissipation  and  mechanical  work  against 
the  large  centrifugal  force  into  the  problem.  The  latter  work  will 
be  much  more  complicated  than  the  present  one.  When  the  angular 
velocity  of  the  system  is  reduced,  the  Coriolis  force  will  gradually 
increase  in  significance:  this  will  results  in  an  asymmetric  boundary 
layer  problem,  and  the  quasi-one-dimensional  conduction  equation  may 
no  longer  be  applicable.  When  the  angular  velocity  is  reduced  still 
further,  the  periodic  effect  of  the  gravitational  field  will  give  rise 
to  an  unsteady  problem  and  the  flow  may  then  oscillate  along  the  fin 
with  each  revolution.  Finally,  when  the  rotating  speed  tends  to  zero, 
it  becomes  a  completely  different  problem  of  an  inclined  fin  in  a  gra¬ 


vitational  field. 
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APPENDIX  A 
NORMALIZATION 

Under  the  assumptions  made  in  section  2.1-1,  the  governing 
equations  of  this  rotating  fin  problem  are: 
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Introducing  the  following  normalized  variables  into  above  equations 
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where  subscript  c  indicates  the  corresponding  characteristic  values  which 
will  be  determined  later.  The  normalized  equations  are  therefore 
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If  the  system  rotates  at  a  constant  high  speed,  the  periodically 

changing  effect  of  gravitational  force  can  be  suppressed  in  relation  to 

the  centrifugal  force  in  equation  (A. 6).  Hence,  when  the  Froude  number 

q 

(Fr) ,  defined  as  R  co  /g,  is  much  greater  than  unity,  and  the  boundary 
conditions  are  independent  of  time,  the  problem  becomes  steady  after  a 
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long  period  of  time,  i.e.,  t  is  large.  In  order  to  compare  the 

coefficients  of  each  normalized  terms  governing  this  quasi-steady 

•  % 

problem,  six  characteristic  values  have  to  be  obtained.  Ab  initio 
L  and  may  be  chosen  as  X  and  0^  respectively.  The  others 

have  to  be  determined  with  reference  to  the  physical  characteristics 
of  the  system. 

In  equation  (A. 5),  8©c  is  normally  small  for  free  convection 
problems  and  the  right-hand-side  is  therefore  negligible,  especially 
when  t  is  large.  Hence  there  are  only  two  terms  left  and  their 
coefficients  must  be  of  the  same  order  of  magnitude, 
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For  convenience  it  can  be  set 
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Since  the  momentum  boundary  layer  in  free  convection  is  manifest 
through  the  effect  of  viscosity  in  the  presence  of  buoyancy,  the  co¬ 
efficients  of  the  largest  viscosity  term  and  the  centrifugal  force 
term  in  equation  (A. 6)  can  be  set  equal.  Hence 
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Assuming  that  advection  and  conduction  in  the  energy  equation 


(A.  8)  are  equally  important  yields 
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and 


From  (A. 9),  (A. 10)  and  (A. 11),  one  can  have 

Y  =  L/Ra1/4 
c 

1/2 

U  =  k  Ra  '  /L 
c 

V  =  k  Ra1/4/L 
c 


where  the  Rayleigh  number,  Ra,  is  defined  as 


Ra  = 


BR  w2  (T  -  T  )  L3 
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VK 


In  the  lateral  momentum  equation  (A. 7),  the  buoyancy  effect  of 
the  gravitional  force  is  suppressed  by  that  of  the  Coriolis  force  when 
the  rotating  speed  is  high,  hence  the  possible  driving  forces  to  cause 
the  pressure  gradient  are  the  Coriolis  force  and  the  centrifugal  force, 
In  order  to  find  which  of  them  is  more  dominant,  the  magnitudes  of 
these  two  terms  should  be  compared.  Their  rotio  is 
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If  Prandtl  number  is  of  order  one,  the  above  ratio  is  generally 


much  greater  than  unity  as  the  rotating  speed  is  high.  Therefore, 


Coriolis  force,  2u<j>,  gives  rise  to  the  lateral  pressure  gradient 

and  it  is  convenient  to  set 

P  =  ptoY  U  30 
c  c  c  c  , 
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Once  the  characteristic  values  are  established,  the  normalized 
equations  for  this  steady  problem  can  be  rewritten  as 
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«  0[1]  ,  (A. 14) 


(A. 15) 


When  air  is  the  fluid,  Prandtl  number  is  of  order  one  and  the 

-1/4  -1 

conditions  that  Ra  ,  Fr  ,  3(T  -  T  )  and  Os  are  all  much  less  than 

7  L  00 

unity,  the  boundary  layer  equations  are  simplified: 


and 
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(A. 16) 
(A. 17) 

(A. 18) 
(A. 19) 


As  an  example,  for  a  6-inch  aluminum  fin,  12  inches  from  the 
center  of  rotation,  rotating  synchronously  with  air  at  300  rpm,  if 


. 


(T  -  T  )  is  10°  F,  one  will  have  Fr  1  -  0.02,  3(T  -  T  )  -  0.019, 

■£  00  7  ?  X  IT  CO  7 

-4  -1/4  -1  -1/2 

Os  -  2.2  x  10  ,  Ra  '  -  0.011  and  (T  -  T  )  Ek  Ra  -  0.1. 

-£  OO 

In  all  cases,  it  is  now  clear  that  the  nature  of  this  quasi¬ 
steady  problem  is  governed  by  the  simultaneous  equations  (A. 16), 

(A. 17)  and  (A. 19),  while  the  separatable  equation  (A. 18)  can  indicat 
the  pressure  distribution,  if  it  is  required,  after  the  velocity 
and  temperature  profiles  have  been  obtained. 

The  normalized  boundary  conditions  associated  with  these 
simultaneous  governing  equations  are 

x  =  0  :  u  =  0,  <J>  =  0 

y  -  0  (x  >  0)  :  u  -  0  ,  v  =  0  ,  <f>  =  4>o  (x) 

y=00:  u  =  0  ,  =  0 

V7here  c|> 0  (x)  is  the  temperature  distribution  along  the  fin  surface. 


, 


52 


APPENDIX  B 

LOCAL  SIMILARITY  TRANSFORMATION 


The  problem  is  to  solve  three  partial  differential  equations, 
(A. 16),  (A. 17)  and  (A. 19),  simultaneously.  Introducing  a  stream 
function  ^(x,y)  defined  by 


dijj  ,  dip 

U  ”  3?  md  V  =  -  37  ,  (B.l) 

equation  (A. 16)  is  automatically  satisfied  and  the  governing  equations 

reduce  into  two  partial  differential  equations: 
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3  3^  3^ 

~  3  3x  3y 
3y 


3^ 

9y 


d2ip 

3x3y 


+  Pr 


(1  -  vx)  4> 


=  0  , 


(B.2) 


d2(p  _  3(f) 

_  2  3y  3x  +  3x  3y 
3y 


The  corresponding  boundary  conditions  become 


(B.3) 


3^ 


=  ° , 
3y 


4»  =  0  , 

(p  =  4>q  (x)  , 

(p  =  0 

where  i)i  =  0  at  y  =  0  is  reasonable,  since  for  boundary  layer  problems 
without  mass  transfer  through  the  wall,  it  can  always  be  arranged  in 


x  =  °  ;  =  0  > 

y  =  0  (x  >  0)  ;  iJj  =  0  , 

dlP  n 

y  =°d  ;  37  “  0  ’ 


this  way  that  the  stream  function  will  be  zero  at  the  wall. 


- 
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The  ensuing  effort  is  normally  in  finding  a  similarity  solution 
for  suitable  form  of  <j)g  (x)  .  However,  it  is  vital  here  to  emphasize 
that  4> o  (x)  is  an  unknown  temperature  distribution  along  the  fin  sur¬ 
face.  It  is  neither  feasible  nor  desirable  to  assume  that  <J>q  (x)  will 
definitely  fulfil  the  special  power  law  requirement  of  similarity. 
Besides,  whenever  4>q  (0)  is  not  zero,  there  is  a  singularity  at  the 
leading  edge,  i.e., 


lim  <j)g(x)  ^  iira  <j>(0,y)  =  0  . 
xrH3  y+Q 

Hence,  it  is  now  necessary  to  seek  a  transformation,  without  prescribing 
the  interfacial  temperature  4>q(x),  to  remove  this  singularity.  A  trans¬ 
formation  can  be  assumed  in  the  form: 


f  (5,n)  = 


Mx,y) 


$ 


U,n;  *0OO  ’ 


(B.4) 


F(x) 

ri  =  H(x)  I(y)  ,  and  5  =  K(x)  . 

where  $(£,ri)  and  f(^,q)  are  two  dependent  variables  in  the  new  coordi¬ 
nate  system  corresponding  to  <J>(x,y)  and  ip(x,y)  respectively. 

The  present  object  is  to  define  suitable  functions,  F(x),  H(x) ,  K(x) 
and  I(y),  to  facilitate  the  numerical  work.  Substitution  of  the  above 
transformation  into  equations  (B.2)  and  (B.3)  followed  by  rearrangement 
give 
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(B.5) 
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and 


(^H2^2)*"  +  O0HIyy)4>' 
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(4»0HF  I  )f$'  -  (FHcfjn  I  )f'$ 


x  y 


x  y 


-  (FHt0KxIy)(f'|i-.-|f)  (B.6) 

where  the  primes  denote  differentiations  with  respect  to  n ;  x  and  y 
derivatives  are  indicated  with  the  subscripts.  It  should  be  noted 
that  both  equations  have  been  arranged  in  such  a  way  that  any  term 
dependent  on  E,  is  kept  on  the  righ-hand-side  of  the  equation.  This 
is  the  essential  treatment  of  constructing  a  basis  for  "local  simi¬ 
larity"  . 

Since  the  highest  derivative  in  the  differential  equation  must 

not  be  eliminated,  equations  (B.5)  and  (B.6)  can  be  simplified  by 

3  3  2  2 

dividing  by  PrFH  I  and  4>qH  1^  respectively.  Hence  the  coefficients 
of  f  ’  '  ’  and  $ ' '  become  unity  in  the  equations:  thit  is,  the  equations 
become 
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and 
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(B.8) 

The  first  and  the  second  square-bracketed  terms  of  equation  (B.7) 
are  related  to  viscous  and  inertia  forces  while  the  first  and  the 
second  square-bracketed  terms  of  equation  (B.8)  correspond  to  con¬ 
duction  and  advection.  The  third  square-bracketed  term  in  equation 
(B.7)  is  the  buoyancy  effect. 

It  is  important  that  in  determining  the  four  disposable 
functions,  F(x),  H(x),  K(x)  and  I(y),  the  pitsfalls  of  eliminating 
the  highest  order  derivatives  in  the  physically  significant  brackets 
of  equations  (B.7)  and  (B.8)  should  be  avoided  and  the  singularity 
of  the  problem  should  be  removed.  Each  disposable  function  can  be 
determined  either  by  directly  defining  it  according  to  one’s  special 
purpose  for  solving  the  problem,  or  by  equating  a  specific  value  to 
the  parenthesized  function  group  of  any  significant  term  in  equations 
(B.7)  or  (B.8) 

The  normalization  process  implies  that  inertia  force,  viscous 
force,  and  buoyancy  force  terms  are  equally  important  in  equation  (A. 17) 
and  have  to  be  retained  in  this  problem.  Besides,  the  advection  terms 
must  be  as  significant  as  the  conduction  term  in  equation  (A. 19). 

Thus,  two  parenthesized  groups  can  be  chosen  and  set  to  be  two  undefined 
non-zero  constants  and  C 


' 

I 


56 


—  =  C 
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Pv-earrangement  of  (B.9)  yields 
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(B.ll) 


and  hence. 
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Choos 


ing  and  =  0,  I(y)  can  be  kept  in  the  simplest  form, 


i  .  e. 


Ky)  =  y  . 


(B . 12) 


From  equations  (B.10),  (B.ll)  and  (B.12), 

<f>o(x)  =  C.F(x)H3(x)  , 
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and 
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(B • 14) 


where  t  is  the  dummy  variable. 

Since  the  transformation  is  to  remove  the  singularity  which 

arises  when  lim  4>o(x)  ^  0  while  lim  4>(0,y)  =0,  it  is  wished  that 
x-*0  y^O 

this  sigular  point  (0,0)  in  x  -  y  plane  can  be  transformed  into  a 
line  in  the  new  S~h  plane,  so  that  this  sudden  change  at  the  singular 
point  becomes  smooth  change  along  the  line.  This  is  fulfilled  by 
setting 


C  ^  =  C  2  —  1  and  = 


1/3 


(t)dt 


in  H(x).  Hence  F(x)  and  H(x)  are  simplified  to  be 
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(B . 16) 


Substituting  equations  (B.12),  (B.15),  (B.16)  in  (B.7)  and 
(B.8),  and  retaining  those  terms  with  coefficients  which  are  functions 
of  x  on  the  right-hand-side  of  each  equation,  yields 
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and 
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(B . 18) 


Since  these  transformed  equations  are  now  in  5-q  coordinates, 
those  coefficients  contained  by  brackets  in  equations  (B.17)  and 
(B.18)  should  be  either  constants  or  functions  of  5  and  n .  The  last 
term  of  equation  (B.17)  suggests  that 

5  =  K(x)  =  yx  .  (B . 19) 


Hence, 


dK(x)  = 
dx  Y  ’ 

and  for  a  given  function  Gq(C),  if  Gp(?) 


<|>o(x)  then 


d4>o (x)  _  dGp (g) 
dx  dC 


and 


x 

$0  (t)dt 


G0(S)dS 


where  Gp  (O  does  not  exhibit  the  same  functional  form  as  4>p(x);  it 
only  implies  that  it  is  the  same  physical  quantity  but  represented 
in  different  coordinate  systems.  By  substitution,  all  the  bracketed 
coefficients  in  equations  (B.17)  and  (B.18)  become  functions  of  K. 

It  is  noteworthy  that  the  salient  merit  of  definition  (B.19)  not  only 
transforms  the  equations  completely  into  the  (C,n)  domain,  but  also 
saves  a  great  deal  of  computing  time.  This  is  because  y  is  never 
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greater  than  one  and  the  number  of  steps  required  in  the  approxima¬ 
tion  along  ^-direction  is  less  for  smaller  y  The  detail  will  be 
discussed  later. 

Therefore,  the  definitions 
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transform  governing  equations  into  a  new  coordinate  system  of  ^  and  q , 


for  which. 
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The  associated  boundary  conditions  in  n-direction  can  easily 
be  deduced  from  those  of  equations  (B.3)  and  (B.4):  whence 

n=0;  f  =  0  ,  f*  =  0  ,  $  =  1  , 

n  =00  ;  f'  =  0  ,  $  =  0  .  (B.23) 


dG_1/ 3 

Provided  that  ^  is  finite  and  the  derivatives  of 
bounded  when  5  =0,  the  right-hand-sides  of  equations 
(B . 22)  vanish  completely  at  this  leading  edge.  Hence 
of  the  boundary  layer  equations. 


f  and  $  are 
(B . 21)  and 
the  solution 


Pr(f'  *  »  +  $)  +  ff  '  '  -  |(f  ')2  =  0  (B .  24) 

and 

+  f$»  =  0  ,  (B . 25) 

and  the  boundary  conditions  of  (B.23)  give  the  boundary  values  along 
the  line  K  =  0. 

Although  Go (5)  is  the  temperature  profile  for  each  of  the 
coupled  heat  transfer  problems,  it  can  be  conceived  of  as  an  unknown 
but  fixed  input  to  the  convection  system  of  the  fin  problem.  If 
Go(?)  is  "known",  the  coefficients  on  the  right-hand-sides  of  equations 
(B.21)  and  (B.22)  can  be  evaluated  at  each  specific  C-level.  In  the 
study  of  the  boundary  layer  problem,  most  interest  is  placed  in  the 
variations  of  velocity  and  temperature  in  the  lateral  direction.  Hence, 
whenever  the  first  partial  derivatives  with  respect  to  5  are  replaced 
with  suitable  approximations,  equations  (B.21)  and  (B.22)  become  two 
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coupled  ordinary  differential  equations  governing  boundary  layers  at 
that  specific  5-level:  this  is  the  local  similarity  concept. 

Several  favorable  features  can  then  be  seen  from  this  "local 
similarity  transformation": 

(1)  The  definitions  of  the  transformation  (B.20)  can  be 
rewritten  as: 

5  =  yx  , 


n  = 


Gn1/3(Q 


~f  G01/3(S)dS 
^0 


1/4 


*Kx,y)  = 


f-f  G0l/3(S)dS 


3/4 


f(5,n)  , 


<Kx,y)  =  GqCO 

and 

4> o  (x)  =  Gq  (5)  . 

The  singularity  caused  by  <j>o(0)  “  Gq(0)  *  0  at  the  leading 
edge  is  obviously  removed,  after  the  point  of  discontinuity 
(x  =  0,  y  =  0)  has  been  expanded  into  the  line  5  =  0  in  the 
5-q  plane.  Along  this  line,  $(0,n)  changes  from  unity 
(q  =  0)  to  zero  (n  ->-oo  ). 

(2)  As  stated  before,  integration  of  equations  (B.24)  and  (B.25) 
which  are  subjected  to  the  boundary  conditions  (B.23)  will 
give  the  supplemental  boundary  conditions  along  line  5  = 
Hence,  there  is  no  difficulty  in  starting  the  numerical 
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integration  of  the  boundary  layer  equations,  (B.21)  and 
(B.22),  in  the  £-r|  plane. 

(3)  It  is  not  subjected  to  the  difficulty  of  deciding  which 
similarity  solution  should  be  used,  since  no  specific 
similarity  restriction  on  the  wall  surface  temperature 
is  of  concern  in  the  transformation.  This  allows  the 
possibility  of  matching  the  interfacial  temperature  Gq(0 
of  such  a  coupled  heat  transfer  problem  of  conduction  and 
convection. 

(4)  In  a  boundary  layer  problem,  comparatively  less  interest 
is  taken  along  the  longitudinal  direction.  This  trans¬ 
formation  introduced  only  first  derivatives  with  respect 
to  0and  this  is  very  helpful  in  applying  the  local 
similarity  concept. 

(5)  All  the  boundary  conditions  are  properly  transformed. 
Hence,  there  is  no  violation  between  the  original  and  the 
transformed  boundary  values. 
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APPENDIX  C 


QUAS I-ONE-DIMEN S TONAL  PROBLEM  OF  FIN 

For  the  simplicity  in  studying  the  validity  of  the  quasi-one- 
dimensional  fin  equation,  a  constant  cross-section  area  and  constant 
property  circular  fin  will  be  assumed.  The  conduction  heat  transfer 


in  the  fin  is  governed  by 


*8To(?,x)*lir  IR 

3X2  8R 


0 


(C.1) 


and  the  boundary  conditions: 


T0(R,0)  =■  q  ,  3T0ax’~  +  [ToO.D  -  Tj  =  0  , 
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Tq(0,X)  =  finite  , 


9Tn(a,X)  ^  h 


9R 


+  [T0(a,X)  -  T J  =  0 

s 


where  L  and  a  are  respectively  the  length  and  radius  of  the  fin,  while 
and  indicate  the  ambient  temperature  and  the  temperature  at  the 
root  of  the  fin. 


After  introducing  the  normalized  geometric  variables: 

R 


X 

x  =  -  , 


r  = 


and  a  new  temperature  variable 


Tn(X,R)  -  T0 
TO  Vx> r  /  t1  —  T 

J.  ±oo 


the  differential  equation  becomes 


9  4>o 

3x2 


+ 


L  7_  _3 
2  8 


r  (r  1^)  ■ 


0 


(C.2) 


*  In  this  appendix,  X  is  the  longitudinal  displacement  from  fin  root, 
and  R  is  the  radial  distance  from  the  axis  of  the  circular  fin. 


. 
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which  is  subjected  to  the  normalized  conditions: 

Mr.O)  =  1  ,  =  Mr>i)  ( 


(0,x)  =  finite  ,  *0(1,*)  . 


Integrating  equation  (C.2)  with  respect  to  r  from  0  to  1  and 


applying  the  last  two  boundary  conditions  give 

}  2  2 
'  3  L  h  A  ^  N 

- u-v  -  r  dr - —  4>o(l,x)  =  0  . 


0 


3x 


aK 


(C.3) 


When  <f>o  is  independent  of  r,  the  problem  is  (by  definition) 
one  dimensional.  When  <J>q  is  almost  constant  along  the  radial  direction 


on  each  cross  section  of  the  fin,  i.e. 

lit 


3r 


«  0[1]  , 


(C  .4) 


the  approximations 

4>0(l,x)  -  (x)  and 


324>n  (r,x) 

,  2 

3x 


d24>n  (x) 

dx2 


may 


be  used  in  equation  (C.3)  to  produce  the  quasi-one-dimensional 


equation 


d  <f> 


dx 


0.  


2L2h  ,  n 

-W  =  °> 


(C .  5) 


or  in  the  original  form 

dx2 


2h 

aK 


(T0  -  T  )  =  0. 


(C.6) 
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In  other  words,  the  necessary  condition  to  construct  the 
quasi-one  -  dimensional  fin  equation  is 


9^0 (r , x) 

9r 


«  0[1] 


Moreover,  since 


9$o(l,x) 

9r 


ha 

K 


$0 (l,x) 


and  $0(1, x>  =  1  > 

one  can  predict  that 


9$0 (r,x)  _ 

9r 

Therefore,  (C.4)  and  (C.7)  infer  that 

Bi  =  —  «  0[1]  (C . 8) 

s 

is  the  necessary  condition  for  the  utility  of  equations  (C.5)  or 
(C.6).  This  conclusion  was  deduced  numerically  [15]  when  Irey  checked 
the  total  amount  of  heat  flux  calculated  from  the  quasi-one-ditnensional 
fin  equation  with  that  of  two  dimensional  Laplace  equation. 

Irey  concluded  his  results: 

"as  far  as  the  more  significant  variable,  the  total  heat 
transfer,  is  concerned,  the  standard  criteria  of  large 
length  to  diameter  ratio  is  not  only  not  of  primary 
significance,  it  is  almost  irrelevent.  Of  far  more 
significance  is  the  ratio  of  interior  to  exterior  resis¬ 
tance,  the  Biot  number,  and  that  only  for  small  Biot 


ha 


K 


(C.7) 
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number  is  the  one-dimensional  solution  a  satisfactory 
approximation. M 

Hence,  when  Bi  <<  0[1],  the  problem  is  governed  by  equations  (G.6) 
and  the  boundary  conditions  which,  in  dimensional  form  are: 

T0(0)  =  T.  and  ^  ^  [T0(L)  -  Tj  =  0  . 

s 

A  related  question  that  frequently  arises  concerns  the  condition 
at  which  the  second  boundary  condition  can  be  replaced  by 


To  a)  =  T 


or 


dTn(L)  =  n 
dX  * 


(C .  9) 


Either  approximation  of  (C.9)  states  that  the  heat  transfer  rate  at 
the  end  of  the  fin  is  zero.  Apparently,  this  is  acceptable  only  when 
the  total  amount  of  heat  dissipated  at  the  end  of  the  fin  is  compara¬ 
tively  small.  Consider  the  situation  when  the  heat  transfer  coeffi¬ 
cient  is  constant,  the  expression  for  the  total  heat  flux  conducted 

away  from  fin  is  given  by 

L 


Q  =  2rah  [T0(X)  -  T J  dX  +  ra  h  [T0(L) 

Jo 


T  ] 


=  2-rraLh  j  <j>0(x)dx  +  ira  h<f>0(l)  , 

0 

it  can  be  seen  that  the  approximation  holds  when 


Tra  h  4*  n  ( 1) 


2iraLh  |  4>n(x)dx 


<<  0 [ 1 ]  , 


0 


' 


* 
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i.  e. 


- gft.Q  (1) -  «  0[1]  . 

2L  J  4> o  (x)dx 

0 


If  and  x  are  normalized,  this  reduces  to 


2L 


«  0[1] 


In  a  word,  if  Bi  <<  0[1]  and  —  <<  0[1],  the  fin  system  is 
approximately  governed  by  the  equation 
d2In(X)  2h 

,2  ~  aK~  lTo(X)  -  T»1  -  0 

Q.X  S 


and  the  associated  boundary  conditions 


T0  (0)  =  Tr  and  T0  (L)  =  To 


“T-  • 


In  the  study  of  conduction  in  a  6-inch  triangular  aluminum  fin 


with  1/2-inch  base  width,  the  empirical  data  cited  in  reference  [16] 

1/4 

1"l) 


give 


h  =  0 . 29  (—^ 


BTU 


hr-ft2-F° 


for  air 


and 


K  =  100  -  200 
s 


BTU 


hr-f  t-F' 


for  aluminum. 


if  AT  =  100°F,  the  maximum  Biot  number  of  the  problem  is 

Bi  =  -  5 . 4  x  10~4  «  0[1] 

K 

s 


and  the  geometric  ratio  is 


Y  «  0.08  «  0 [ 1] 

Li 


' 


It  is  apparent  that  in  this  case, 
conduction  equation  and  the  approximate 


the  quasi-one-dimensional 
boundary  condition  (C.9) 


are  acceptable. 
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APPENDIX  D 


FORMULATION  OF  THE  LIMITING  CASE  y  =  0 


From  appendix  B,  the  similarity  transformations  are 

5  =  Yx  , 

1/3, 


n  = 


$0  (x) 


X 


3  f  *01/3(S)dS 


'0 


1/4 


r  x 


^(x,y)  = 


4[  +o1/3(s)ds 


3/4 


f(5,n)  , 


and  <Kx,y)  =  4>o(x)  (K5>p) 

As  y  ->  0 ,  the  variable  5  ■>  0j  a  small  modification  of  this  transform 
mation  for  this  case' is  to  set  £  =  x,  and  then  the  boundary  layer 


equations  become 

Pr  (f,,T  +  $)  +  ff"  -  ^  H1(4»0)  ( f  * ) 2  +  H2((])o)(f,|f 


-  f’’-^) 

35/ 


+  Pr  y5$ 


and 

$' '  +  f$f 


=  3H  (<j>0)f  '*  +  H2(4>0)  (f 'll  +  *'ff) 


(D.l) 


(D.2) 


where  and  H2  are  defined  as 


H-.  (4>0)  =  - 


4  dtn~1/300  f  4,01/3(s)ds 
3  dx  1 


(D.3) 


and 


H0  (<|>o 


)  =  3  4>0-1/3(x)["  to1/30)ds 


(D.4) 


, 
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It  is  obvious  that  as  y  ->  0,  the  last  terra  in  equation  (D.l) 
becomes  zero  while  H^(4>0)  and  maY  not  vanish.  Since  £  =  x. 

Go  (?)  =  ^(x)*  i.e.  not  only  numerically  equal  but  also  of  the  same 
form.  Thus  as  y  ->-  0,  equations  (D.l)  and  (D.2)  become 


Pr  (f"  +  *)  +  ff"  (f’)2  =  H  (G0)(f') 

+  h2(g0) 


.n2 


.JT 

3? 


fMl 

3?J 


(D.5) 


and 


+  f$’  =  3H1(G0)ff$  +  H2(G0)Ff,||  -  $’||]  , 

the  transformed  boundary  conditions  are 


(D.6) 


q  =  0;  f  =  f'  =$-1  =  0 

n  =  oo  ;  f '  =  $  =  0, 


and  the  solution  of  this  boundary  value  problem  at  the  leading  edge 
will  give  the  "initial"  condition  at  ^  =  0.  Similarly,  the  conduction 
in  solid  is  governed  by 


d_  (  F  dGg(Q\ 

d?  d?  / 


+  C<D'(?,0) 


jf  G01/3(S)dS 
d? 


3/4 


G0 (?)  =  0 


(D.7) 


wi  th 


dG  n  (  ?  )' 


d? 


=  0 


and 


G0 (1)  =  1 


?=0 


Equations  (D.5)  and  (D.6)  indicate  that  in  spite  of  the  necessity 
of  satisfying  the  above  conduction  equations  and  their  associated 


- 
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boundary  conditions,  the  similarity  solutions,  if  they  exist,  should 
also  meet  the  following  two  conditions: 

(1)  function  H^(Gq)  should  be  independent  of  and 

(2)  either  function  H^CGq)  is  zero  everywhere  or  §  and  f  are 
function  of  n  only  (with  H^Gq)  bounded). 

In  the  first  case, 

Hl(Go)  =  ~  ~f  ~"°  d£  (0  f  G01/3(S)dS  =  A1  ,  say,  (D.8) 

J0 

and 

H2(Go)  =  -f  G0"1/3(O  I  G01/3(S)dS  =  0  ,  (D . 9) 

•^0 

when  condition  (D.8)  is  integrated  with  respect  to  £,  one  can  easily 

4 

verify  that  =  —  with  the  help  of  condition  (D.9).  Whereas,  no 
wall  temperature  Gq(5)  can  satisfy  condition  (D.9)  and  hence  similarity 
solution  is  impossible  in  this  case. 

In  the  second  case, 

-1/3 

VGo)  =  -  |  dG°  d;  I  G01/3(S)dS  =  A2  ,  say,  (D.9a) 

J0 

and 

$  =  $(ri)  and  f  =  f(n).  (D.10) 

Condition  (D.10)  implies  that  those  bracketed  terms  in  equation  (D.5) 
and  (D . 6)  vanish  automatically.  Besides,  in  conduction  equation  (D.7), 
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$’(£»0)  ls  no  longer  a  function  of  £,  i.e.  (£,0)  =  $'(0).  Hence 
the  possibility  of  similarity  solutions  relies  on  the  satisfaction 
of  (D.8)  and  the  conduction  problem. 

It  can  be  shown  by  substitution  of  the  general  power  law  surface 
temperature  assumption,  Gq(£)  =  (a  +  b£)n,  into  the  conduction  part 
that  only  Gq(£)  =  £  may  present  similarity  solution.  If  this  power 
law  assumption  is  made.,  then 

VG»>  ’  J^rf-3) 


and  equation  (D.7)  becomes 


2  n-1 

n  £ 


15n  -  3 
12 


0  . 


(D .  11) 


The  fact  that  this  equation  must  be  satisfied  at  every  value  of  £ 
between  zero  and  unity  requires  that  the  power  of  £  in  both  terms 
must  be  equal  and  then  the  sum  of  the  two  coefficients  must  be  equated 
zero.  Unfortunately,  the  first  requirement  yields  n  =  -3  and  with 
this  value,  not  only  can  the  second  requirement  not  be  satisfied  but 
also  the  boundary  condition  at  the  tip  of  the  fin  is  inadmissible. 
However,  there  is  one  special  case  for  a  similarity  solution.  When 
n  =  0,  i.e.  an  isothermal  surface,  and  C  =  0,  equation  (D.ll)  is  always 
true  and  H  (Gq)  =0.  No  conflict  exists  between  conditions  n  =  0  and 
C  =  0,  on  the  contrary,  the  former  is  a  consequence  of  the  latter. 


!  I? 


In  chapter  II,  C  is  defined  as 


C  =  k  h  1/4 
L  W  K  R 


s 

As  far  as  a  triangular  fin  is  concerned,  the  slenderness  pj-  will  not 

w 

likely  be  zero.  The  situations  for  C  =  0  are  therefore: 

(1)  K  =  0,  i.e.  insulated  fin  surface,  and  obviously  it  will 
give  rise  to  an  uniformly  heated  body  (n  =  0) ,  where 


T  =  T  . 
r 

(2)  K  =<X>,  although  physically  a  solid  with  infinite  thermal 

s 

conductivity  does  not  exist,  it  still  can  be  conceived 
that  such  a  perfect  conductor  will  be  an  isothermal  body 
under  steady  state  condition. 

(3)  Ra  =  0  due  to  0  =  0.  Apparently  when  T^_  =  T^,  the  entire 

fin  is  of  the  same  temperature  T  . 

In  fact,  the  existence  of  an  isothermal  fin  for  C  =  0  can  be 
used  when  y  k  0.  As  C  =  0,  the  conduction  equation  (D.7)  reduces  to 


After  integration  and  using  boundary  condition 


it  yields 


d G 
d£ 


A-  =  0 


one  more  integration  and  the  isothermal  profile  is  obtained  following 


n 


, 


the  application  of  G0(y)  =  1.  Hence,  when  C  =  0 


G0(O  =  1. 


It  is  evident  that  this  solution  is,  in  a  sense,  trivial.  No 
matter  whether  C  =  0  because  of  situation  (1)  or  (3),  no  free  con¬ 
vection  can  take  place.  However,  it  can  be  a  good  approximation  to 
the  solution  for  convective  fin  with  very  high  conductivity  in  the 
solid.  In  addition  to  this  practical  usage,  the  behavior  of  the 
isothermal  surface  can  be  used  as  a  way  of  checking  the  "local  simi¬ 
larity  transformation"  of  this  study.  Mien  Go  =  1  throughout  the 
fin,  the  transformation  becomes 


5  = 


x 


x 


y 


1/4 


and 


These  are  the  same  forms  as  were  developed  for  the  classical  isothermal 
plate  problem  in  references  [1,3, 4, 5]. 
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APPENDIX  E 


NUMERICAL  TECHNIQUE 


It  has  already  been  shown  that  this  rotating  fin  problem  is 
governed  by  the  boundary  layer  equations; 

Pr(f,,f  +  $)  +  ff”  -  4(f')2  =  F  (G0,O(f')2 


+  F2(G0,5)  ff'll’-  f"||)  +  Pr5$  ,  (E.l) 


+  f$'  =  3F1(G0,C)f'$  +  F2(G0,C)  (f|§  -  $'||) 
and  the  quasi-one-dimensional  equation  of  conduction. 


(E.2) 


dsUdfj  +  ^3^0  >  5)Go  -  0  , 

with  the  associated  boundary  conditions; 


and 


provided  that 


dG 


n  =  0 

q  =00 

5  =  o 

5  =  Y 

-1/3 


f=f'=$-l=0, 

f'  =  $  =  0  , 
rdGa  =  o 

Si  5  * 

G0  =  1  » 


(E.3) 


(E.4) 


(E.5) 


d? 


is  finite  and  the  derivatives  of  f  and  $  are 


bounded  when  5=0.  Functions  F^  F2  and  F3  are  defined  as 


Ft  (G0,  O 


4  „  -4/ 3,  d GqO)_I  g  1/3 


T  G0 


(0 


d5  J  ~° 

o 


(S)dS 


F2(G0,O  =  j  G0"1/3(O J  G01/3(S)dS  , 


0 
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and 


f3(g0,O 


=  C$'(5,0) 


d_ 

d? 


1/3 

G0  /J(S)dS 


3/4 


The  complexity  of  these  three  coupled  governing  equations 
reveals  that  the  analytic  solution  can  not  be  procured  by  means  of 
established  mathematical  techniques.  A  numerical  method  is  hence 
considered  as  the  most  suitable  scheme  in  solving  this  problem. 

Apparently,  the  conduction  equation  is  coupled  with  the 
boundary  layer  equations  through  F^,  F^  and  F^  which  are  functions 
of  the  interfacial  temperature  distribution  Gq(C).  Therefore,  the 
matching  problem  between  conduction  and  convection  systems  can  be 
solved  with  a  proper  iterative  routine  searching  for  the  true 
interfacial  temperature  Gq(6). 

From  the  standpoint  of  free  convection  alone,  the  surface 
temperature  Gq(C)  can  be  regarded  as  the  input  of  the  problem.  With 
a  reasonable  guessed  profile  G0(Q  increasing  from  a  positive  value 
at  ?  =  0  to  unity  at  5  =  Y,  the  strongly  coupled  boundary  layer 
equations  can  be  solved  using  the  concept  of  'local  similarity". 

As  stated  in  appendix  B,  the  problem  of  starting  the  integration 
is  overcome  by  choosing  the  first  5~ level  at  the  leading  edge  where 
both  functions  F^^  and  F2  vanish  (as  5  =  0),  and  equations  (E.l)  and 
(E.2)  become 


Pr(f1 


+  v  +  fi£i"  -!(v)2  =  0  ■ 


i » » 


- 

, 
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and 


$l' '  +  f  $  *  =  0 


y 


where  the  subscript  denotes  the  ordinal  number  of  £-level. 

The  region  is  then  divided  into  several  equally  spaced  ^-levels 
in  order  to  have  an  easier  estimation  of  the  truncation  errors  of  the 
approximations  in  ^-direction.  Besides,  it  can  decrease  the  computing 
time  and  roundoff  errors  on  evaluating  the  non-similar  terms  on  the 
right-hand-sides  of  equations  (E.l)  and  (E.2).  Before  the  numerical 
work  is  applied  to  the  successive  ^-levels,  some  numerical  approximations 
of  equally  spaced  interval  with  their  associated  truncation  errors  should 
first  be  introduced:  if  I<  is  a  function  of  the  independent  variable  S, 
and  AS  is  the  step  size,  then 

(1)  two-point  backward  finite  difference  formula 


3K 


n 


3S 


=  Yq  (K  -  K  )  +  o[|^|  , 

AS  n  n-1  L2  J 


(2)  three-point  backward  finite  difference  formula 


3K 


n 


3S 


2 

(3K  -  4K  .  +  K  0)  +  0[AU 

2AS  n  n-1  n-2  L  3  J  * 


(3)  central-difference  method 


3K 


n 


3S  2AS 


<Kn+i  -  K„-i>  +  °[if  ] 


(4)  two-point  trapezoidal  rule 
S 


f  f  (Kn  +  Kn_p  +  0[f|  ]  , 


n-1 


I 
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(5)  Simpson’s  rule  for  even  number  of  intervals,  i.e.  n  =  odd, 
S 


f  n  As  V* 

K(t)  dt  =  —  )  (K  +  4K  +  K  .)  +  0 
J  3  /  ,  m  m-1  m-2 

0  m=3 


AS 


180 


At  the  second  5-level,  two-point  haclcward  finite  difference  is 
employed  in  representing  the  ^-derivatives  of  f,  f’,  and  $,  while 


dG q.  , 


d5 


is  approximated  with  central-difference  method.  The  integral  of 


1/3 

Go  is  substituted  by  the  two-point  trapezoidal  rule.  The  ordinary 
differential  equations  to  be  integrated  with  respect  to  q  at  this 
level  become 

Pr(f2"'  +  i2)  +  £2f2"  =?G2~4/3<G3  -  V(G21/3  +  Gl1/3)(f2')2 


1  +  (G1/G2) 


1/31 


(f2,}  "  VV  “  f2"f2  +  f2,,fl 


-I-  Pr?2$2  +  0[A5  ] 


and 


>2"  +  f2*2'  =  3  g2  4/3(g3  -  gx)(g21/3  +  G11/3)fCi>. 


2  2 


+  ! 


1  +  (G1/G2) 


1/3 


(f2,$2  "  f2,$l  "  Vf2  +  W  +  °[AC2] 


The  ^-derivatives  of  f,  f'  and  $  in  the  successive  levels  are 
approximated  using  the  three— point  backward  finite  difference.  fhe 
central  difference  is  still  applied  to  ^  throughout  the  remaining 
^-levels  but  the  last  one  where  is  substituted  by  three-point 

backward  finite  difference.  At  each  odd  5-level  (with  an  even 
number  of  intervals),  the  integral  of  Gq1/3  is  evaluated  with  Simpson's 
rule,  while  at  the  even  5_leveis>  integral  over  the  last  interval 


I 
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of  £  can  be  replaced  by  the  two-point  trapezoidal  rule.  Therefore, 

it  is  easy  to  show  that  if  n  is  odd  and  <  £  <  Y, 

Z  n 


Pr(f  ' ’ '  +  $  )  +  f  f  ' '  -  |(f  ')2 
n  n  n  n  3  n 


h  Gn'V3(Gn+l  -  °n-d 


n 


i-  n 


I 

_m=3 


(G  1/3  +  4G  1/3  +  G  1/3) 
m  m-1  m-2 


(fn'>2 


+  |  G  -1/3 
9  n 


Y  (G  1/3  +  4G  1/3  +  G  1/3) 
/  ,  m  m-1  m-2 


m=3 


[V 


(3f  '  -  4f  '  +  f  ') 
n  n-1  n-2 


r\ 

-  f  "(3f  -  4f  .  +  f  0)1  +  Pr£  $  +  0[A5  ] 

n  n  n-1  n-2  J  n  n 


and 

$  ?  1  +  f  $  ’  =  G  "1/3(G  ,  ,  -  G  ) 
n  n  n  9  n  n+1  n-1 


r  n 


I 

m=3 


(G  1/3+4G  1/3+G  1/3) 

m  m-1  m-2 


f  ’$ 
n  n 


r~  n 


+  |g-1/3 
9  n 


V  (G  1/3  +  4G  1/3  +  G  1/3) 

/  ,  m  m-1  m-2 

m=3 


r f  *  (3$  -4$  1  +  $  9) 

L  n  n  n-1  n-2 


-$  ?  (3f  -  4f  -.  +  f  0)1  +  0[A£3] 

n  n  n-1  n-2  J 


If  n  is  even. 


Pr(fn’"+»n)  +  fnfn”  -f(V)2 


2_  -4/3,  _r  , 

(Gn+l  <A-1) 


n-1 


ty  (G  1/3H-4G  .1/3«m  21/3)+  |(Gn1/3+Gn-l1/3) 

3  /  i  v  m  m-1  m-2  z.  n  n  r 


m=3 


2 

(f  ') 

n 


+  |g~1/3 

3  n 


IVCO  1/3+  G  ,1/3+Gm  ,1/3 

3Z_>  m  m-1  m-z 

m=3 


1  ,  1/3  1/3, 

+  2  (Gn  Wn-1 


[ 


f  '  (3f  ’ -4f  /+f  ’) 

•—  v  «  pi— 1  n— 2 


n  n  n- 


-  £n”  (3fn-4fn-l+£n-2)]  +  Pr5A  +  °lA52] 


' 
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and 

$  "  +  f  $  ’  =  G  “4/3(G  l  -  G  J 
n  n  n  3  n  v  n+1  n-1 


47(0  1/3  +4G  3/3  +  G  ,1/3) 

3Z_j  m  m-1  ra-2 


+  y(G  1/3  +  G  1/3) 
2  n  n-1 


+  7Gn1/3  +  Gn  l1/3) 


m=3 
r  n-1 


f  ’$  +  4  G  1^3 

n  n  3  n 


4V  (G  1/3+4G  1/3+G  1/3) 

3/  ,  m  m-1  m-2 

_  m=3 


f  ' (3$  -  4$  .  +  $  „) 

n  n  n-1  n-2 


-  V(3fn~  4£n-l+  fn-2>]  +  0^2> 


At  the  last  5-level,  5  -  Y  and  G  =  1,  if  n  is  odd, 

n 


Pr (f  ’  '  '  +  $  )  +  f  f  •  '-  4  (f  ')2 
n  n  n  n  3  n 


(3-4 G  -  +  G  ) 
27  n-1  n-2 


r-  n 


E 


(G  1^3  +  4G  1/3  +  G  1/3) 
m  m-1  m-2 


+ 


n 


I 

_m=3 


_m=3 

(G  1/3  +  4G  1/3  +  G  1/3) 
m  m-1  m-2 


(f  o2 

n 


\f  1 (3f  ’  -  4f  '  +  f  ’) 
[  n  n  n-1  n-2 


-  f  "(3f  -  4f  ,  +  f  0)1  +  Pry3>  +  0[A53] 

n  n  n-1  n-2  J  n 


and 


$  f  1  +  f  $  '  =  £(3  -  4G  -i  +  G  ) 

n  n  n  9  n-1  n-2 


r  n 


I 

m=3 


(G  1/3+4G  1/3+g  ?1/3) 

m  m-1  m-2 


f 

n  n 


+ 


n 


E 

,m=3 


(g1/3+4g  1/3  +  l/3) 

v  m  m-1  m-2 


[f  *(3$  -  4$  .  +  <Sf  0) 

L  n  n  n-1  n-2 


-  ft  ' (3f  -  4f  -  +  f  9)1  +  0[A53]  , 

n  n  n-1  n-2  J 
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and  if  n  is  even. 


Pr (f  '  r  f  +  $  )  +  f  f  1  ’  +  —  (f  M2 
n  n  n  n  3  v  n  ' 


g  (3  -  4G  +  G  ) 
y  n-1  n-2 

r  n-1 


+  1 


~  n-1 

£ 

L  m=3 

1/3  ,  „  1/3,  .  1  1/3 


(G  1/3  +  4G  1/3  +  G  01/3)  +  ^G  1/3 
m  m-1  m-2  J  2n-l 


(V>2 


£  (G™1/3  +  4vr  ■ +  e^-."  a  +  * 


m-2 


2  n-1 


_  m=3 


f  T  (3f  '-  4f  /+f  ') 

.  n  n  n-1  n-2 


“  £n"(3£n  -  4£n-l  +  W]  +  Pr^n  + 


and 


$  "+f  $  ’=  |(3-4G  +G  „) 
n  n  n  3  n-1  n-2 


(G  1/3+4G  1/3+G  1/3)+  h  1/3 

3/  ,  m  m-1  m-2  2  n-1 

_  m=3 


f 

n  n 


yy (g  i/3+4g  i/3+g  j-/'j)+^g 

3/  ,  m  m-1  m-2  2  n-1 

„m=3 


1/3, ,  1„  1/3 


Ff  *  (3$  -4$  _+$  _) 

L  n  n  n-1  n-2 


’  (3f  -4f  -+f  _)1  +  0[A?2]  . 

n  n  n-1  n-2  J 


The  coupled  partial  differential  equations  of  the  boundary  layer 

problem  have  thus  been  replaced  by  a  finite  set  of  coupled  ordinary 

2 

differential  equations  with  an  error  of  0[A£  ].  Boundary  condition 
(E.4)  indicates  that  it  is  a  boundary  value  problem  at  each  C-level. 
After  the  values  of  f1'  and  at  n  =  0  have  been  guessed,  the 
integration  with  respect  to  n  can  be  carried  out  with  a  fourth-order 
Runge-Kutta  method.  A  scheme  has  to  be  set  up  to  modify  the  initial 
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guesses  of  f,T  and  $  until  the  boundary  conditions  f'  =  0  and  $  =  0 
outside  the  boundary  layer  are  satisfied. 

A  few  trial  runs  have  shown  that  the  numerical  integration  is 
extremely  sensitive  to  the  initial  guesses  of  f"  and  $'  at  n  =  0. 

The  numerical  routine  is  unstable  if  the  guessed  f''(£,0)  and 
$'(£>0)  are  not  very  close  to  their  correct  values.  The  iterative 
technique  of  Nachtcheim  and  Swigert  [14]  with  a  least-squares 
convergence  criterion  on  satisfying  the  asymptotic  boundary  conditions 
has  been  adopted.  Preliminary  tests  indicates  that  when  the  asymptotic 
boundary  conditions  are  satisfied  at  n  =  2  with  a  condition  on  mean 
square  error  being  less  than  10  the  numerical  integration  is  always 
stable  and  convergence  of  fr,(£,0)  and  $'(5,0)  can  be  assured  with 
their  initial  guesses  being  1.0  and  -1.0  respectively.  Using  the 
converged  f,T(£,0)  and  $’(5,0)  as  the  initial  guesses,  the  integration 
(following  the  same  scheme)  can  subsequently  be  carried  out  to  q  =  4, 

6,  8  and  then  10  without  any  special  difficulty  in  either  the  stability 
or  the  convergence  of  the  numerical  program.  The  only  disadvantage 
of  this  technique  is  the  time  of  computation,  even  with  a  IBM  360/67 
high  speed  computer.  For  saving  computing  time,  large  step  sizes  for 
£  and  q  with  an  acceptable  accuracy  are  required.  The  truncation  error 
of  the  coupled  ordinary  differential  equations  is  shown  to  be  of  the 
order  of  A£2.  In  this  study  the  step  size  A£  is  taken  to  be  0.02  so 
that  the  truncation  error  due  to  the  "localization"  of  the  problem  is 


/ 
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less  than  one  percent.  Tests  were  conducted  with  different  step  sizes 
of  p,  and  the  converged  fT,(5>0)  and  $’(£,0)  were  compared  separately 
at  two  ^“levels.  The  comparison  suggested  that  a  step  size  of  0.04 
for  the  integration  in  p-direction  can  give  an  accuracy  to  these 
eigenvalues  within  an  order  of  one  percent. 

This  numerical  solution  of  the  boundary  layer  equations  provides 
the  values  of  (5,0)  at  each  5~level  for  the  initial  guessed  profile 
of  Gq(5).  These  values  of  $’(5,0)  and  Go  (5)  are  substituted  into  the 
convective  coefficient  F^  along  the  fin  surface.  By  means  of  a 
fourth-order  Runge-Kutta  method,  equation  (E.3)  can  be  integrated 
from  the  root  with  the  boundary  condition  Gq (y)  =  1  and  a  guessed  value 


of  5 


dGq 


dG 

.  to  give  — z-Q-  and  G0(5)  at  each  level,  except  the  leading 

?  '  5=y  ^ 


gives  an 


^  *  /dG  \ 

edge.  When  5  =  0,  the  numerical  integration  of  (‘jp’0') 

\d^  /^=o 

infinite  value  of  Gn(0)  since  the  integrand  itself,  obtained  from 


.dG 


R 


’d5 


5=0/ 

re-written  as 


,  is  infinite.  However,  equation  (E.3)  can  be 


5=0 


5 


d2Ga  +  dG_Q_ 


d5' 


d5 


+  FqG0 


0  . 


(E .  6) 


Integration  of  this  equation  simultaneously  with  equation  (E.3)  using 
the  Runge-Kutta  method  can  respectively  give  the  values  of  G0(5)  and 
^Gq.  at  each  level.  A  shooting  method  is  then  employed  to  modify  the 

d5 

guessed  value  of  5f>)  until  the  other  boundary  condition. 


’<15  ) 


5=0 


5=y 

=  0  is  satisfied  to  five  decimal  places 


*• 

•f 


Logically,  this  newly  obtained  Gq (Q  and  the  former  GqCO 
should  give  a  better  guessed  input  to  the  boundary  layer  system  and 
the  whole  procedure  is  then  iterated  until  Gq (0  converges  to  the 
true  temperature  distribution.  However,  practical  experience  points 
out  that  the  computing  time  in  solving  the  boundary  layer  problem 
is  much  longer  than  the  time  on  integrating  the  fin  equations. 

Hence,  it  is  better  to  refine  G q(£)  through  the  conduction  equations 
as  much  as  possible.  Mien  an  improved  Go  (0  ,  obtained  from  halving 
the  lately  obtained  solution  and  the  previous  values  of  Go  (0  ,  is 
substituted  into  equations  (E.3)  and  (E.6)  to  evaluate  a  new  set  of 
with  the  same  values  of  $'(00),  an  iterative  loop  in  integrating 
these  two  equations  is  formed.  The  finally  converged  Gq(0  in  this 
conduction  loop  is  then  used  as  the  improved  input  of  the  boundary 
layer  problem  to  generate  another  set  values  for  $'(5,0).  Such 
iteration  between  the  conduction  loop  and  the  boundary  layer  problem 
continues  until  Go (0  converges  at  every  5— level  with  an  accuracy  of 
ten  decimal  places. 

In  order  to  reduce  the  round-off  error,  double  precision  was 
used  throughout  the  calculation.  This  requires  a  computing  for  a 
complete  solution  of  about  30  minutes  using  an  IBM  360/67  machine. 
Followed  is  the  computer  program  used  in  this  work: 


. 

' 


ooooo  ooooooo 
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PROGRAM  FOR  SOLVING  COUPLED  HEAT  TRANSFER  PROBLEM 
OF  CONDUCTION  AND  LAMINAR  FREE  CONVECTION  OF  A 
HEATED  TRIANGULAR  FIN  IN  A  CENTRIFUGAL  FORCE  FIELD 


EXTERNAL  V  I ,T I, V 1 1 ,T 1 1 , VI  1 1  ,TI 1 1 
REALMS  GXJ(IOO) ,GXK (1 00) , GX J1 3 1100 J , GXK 13 1100 ) , 
1HEAT ( 100) ,HEAF(  100) ,DIFSUM(  100)  ,GD011 ,GD021 , 

2GD031 ,GAMA43, HI  I,F2I 

REAL* 8  R13,R23,R14,R53,R43, SS , S T , SL , E  I  ( 100 ) , 

IS  I GMA ( 100 ) ,GX 1 13( 100)  ,DIF ( 1 00 ) , RAT{ 100 > ,GDI  (100), 
2GDI 1( 100) , GDI I< 100  I, GDI  1 1(100), GAMA,CN,CL ,GDK 11 , 
3GDK 12  *GDHH  ,  GDHH1 , GDK2 1 »GDK2  2 , GDK31 , GDK32 , GDK41 , 
4GDK42,RESI ,GDCOR,DI ,DII ,GX, GXO( 100) ,GA  MA3 
REAL* 8  SF2( 251) ,SFl(25l ), SF (251), SH 1(251),  SH( 251) , 
1SFGMI 1(251) , SF 1  MI 1(251) ,SHOMI 1(251)  , SH 1MI 1  (  251  ) , 
2SFQMI2I 251 ) ,SF1MI2( 251) , SHO  MI  2 ( 25 1 ) , SH 1MI  2  (  25 1 )  , 
3GXI (100) ,XI ( 100 ) ,AC,BC,CC,DC, EC,FC, GC, HC,PR 
COMMON  SF2,SF1, SF,SH1  ,SH,SFOMI 1,SF1MI 1 ,SH0MI1, 

1  SHI  MU,  SF0MI2  ,SF1MI2,  SHOM  1 2  ,  SH l MI  2  ,  PR  ,  AC  ,  8C  ,CC  , 
2DC,EC,FC,GC,HC,GXI, XI , IX,  I , J ,K0 , KOA , KG B ,KOC , 

3KCD , KOE 

READ  (5,100)  K0,K0A,K08,K0C ,KOD,KOE 

100  FORMAT ( »  ' ,613) 

I  J=  1 

I  L=1 

R13= 1*0 DC/3. ODO 
R 14=0  *2  500 
R23=2 .0 D0/3.0 DO 
R43=l .0D0+R13 
R53=l .0D0+R23 
LGOP= 1 
I TERA= 1 

READ  PARAMETER  C  ( CN )  » PRANDTL  NUMBER  (PR), NUMBER 
OF  DISCRITIZED  XI-LEVELS  ( NM) , AND  THE  GUESSED 
EIGANVALUES, F2 I  AND  H1I,F0R  EVERY  XI-LEVEL 

READ  (5,101)  CN 

101  FORMAT  ('  *,023.16) 

READ  (5,102)  PR,F2  I, HI  I , NM 

102  FORMAT  ('  *,3023.16,  14) 

SIGMAQ  )=Q.0 

D  IF ( 1  )  =  C.0D0 
R  A  T ( 1 )=0 • ODO 
NCDD=NM/2*2-l 


M3  Jh  0  9  H2MAXT  T  A  H  03  J9U03  DHJV  02  30 

1  .  : ' !  ’ 

0J  13  3  0*0r  J  A0U3 1 TV  3  A  A\  H  3  PAJOOHAIflT  03T  A'!  H 


1  T,  :  I  *  0,  I  !  1  ,  T  V,  .  *  t  IV  li  ,  'T).l 

,  (  00  (  )  £IXX3,  (  0  M  )f  ILX ,  (  n  *X0,  (COi  )  1X0  r*.!/;  3  q 
,  ■  ■ ,  ■  l  •  ,1  )i  )  •■■  ■  .  (  .  G  :  )  ' ,  ( 

IS  3,]  IH,0  AAf^  A0#  i  ”  on  os 

«(  M3,  ,  »  ;  ■  .  i  B*Ji 

.  C  00  1  )  if  Of  (  001  )TAS  ,  (Of  I }  310,  (  '01  )E  l  I  XO,  <001  )  A’  0  !  ?A 

- 

,  !Axoo,sexo;  ,  n  >oo*ssvc  ,isxoo,p  <oo,c  , si  *'  £ 

.  -  ,  • 

,  <  i  as  )m ? ,  { i  ? s ) i  f  a ,  <  i  as )  ■  ,  <  i? :  > ;  '  v ,  (  i  as  >  s  ia 

.  I  aS  :  !  I  ,  f  1  ■•■:.  >  !  7  ,  (  !  :  :  '  :  V  ’  >r  '  *  ■  “ 

,  <i?s)  u'iha,<fe  »s  iM  rr-  a,  u  «*s  >si  ioa ,  1 1  isimciss 
h  ,  ■’  ,  •  ,:.'i ,  jcuoo,  v  ,  )oii  i  x,  coo  :  >  :>  ■  • 

,  I1M0H2,  l  !:.  I  ?,  i  1^032  ,H2*  Ih2, 32 ,1^  ,  S  "  IDO 
»  30, 0  , 0 A . fl9,SIMIH2,S!M0H2 *  S I  •  I3> ,  v  032, i I ' '  iH2i 
,  j.  x,  hi  a  ,  a  ■  • ,  r  ,  l  ,  i ,  x  i ,  i  x ,  i  x.  ,oi  ,  or ,  0': ,  0  ,  f 

/  ,  •  )>F 

BOX  ,  00  ,  3VX »  iO>  ,  AC-  ,  OX  <'01,r) 

d!o,'  '  )  T A ' ' q 0 3  COI 

I  =L  I 

I -  J  T 

o<  o.£\:*o( .  i  ri  p 
00  r.0=^J 
or  :.£\  }r  .  =  cr 

FI  +■  '  ••  f  f  ’ 

£Sfl+000.1=£?« 

38MUH  JTOHA«3,(HO>  3  9  3T3MA9A9  0A3fl 

?  •  '  !0T  '  ‘  '  )  a  •  J  -  I  >  :  SI  •  I  V  •:  1  - 

V  -IX  Y*  ]\  3  "03.  li  :  A  ]  S;  ,  UJAVMA01  3 

HO  (IOI,a)  GA3fl 
(dl.fSQ,*  *)  7  A  •  0  1 01 

MH,  1  If  ,1  S3,«q  (SOI, a)  0A3P 

(  I  ,di .(  sr  1 1  ♦  »)  tampo-o  soj 
0.0  =  <])A’  01 2 
0  Q  0 »  0  =  (I  )  ^  J  C 
.(  t  ) TAJ 
I-S^SVMl'-nOOH 
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I V=NGDD+1 
MN=NM-1 
C 

C  READ  THE  INITIALLY  GUESSED  INTERFACIAL  TEMPERATURE 

C  DISTRIBUTION, GXI,  AT  EACH  XI-LEVEL  AND  CALCULATE 

C  THE  COEFFICIENTS  REQUIRED  IN  THE  INTEGRATION  OF 

C  BOUNDARY  LAYER  EQUATIONS 

C 

READ  (5,103)  GXim,XI(l) 

103  FORMAT  ( »  *  ,2023.16) 

DC  I  I Z=2 , NM 

READ  (5,103)  GX I ( I Z ) , X  I (I l ) 

DIF(IZ)=XI( IZ J-XI ( 12-1) 

RAT (I Z )  =  D  I  F  ( IZ3/XICIZ) 

1  CONTINUE 
GAM A=X I ( NM ) 

GAMA3=GAMA**3*R43 
GAMA4  3=R43/GA  MA 

2  I  K=  1 

GXO ( I ) =GX I (1) 

0X113(1  )=GX I (  1) **R13 
DC  3  IZ=2 » NM 
GXO ( I Z } =GX I ( I Z ) 

GXJ ( IZ— 1 >=GXI ( IZ-1)+(GXI( IZ )— GX I ( IZ-1 ) )/3.0D0 
GXK( I Z — 1 )=GXJ ( IZ-1 )  +  (GXI( IZ  )-GXI(  IZ-1) )/3.0D0 
GXI 13 ( I Z ) =GX I « IZ)**R13 
GXJ 13 ( IZ-1 )=GXJ (IZ-1)**R13 
GXK13 (IZ-1 )=GXK< IZ-1 )**R13 

3  CONTINUE 

DC  4  1=1,251 
SF1MIK  I  )  =  0 .000 
SFOMI 1( I ) =0.000 
SH1MIH  I  )  =  0 .000 
SHOMI 1( I ) =C. ODO 

4  CONTINUE 
SS=0. ODO 

c 

c  START  SOLVING  THE  BOUNDARY  LAYER  PROBLEM  WITH 

C  SUBROUTINE  "LOSIM"  TO  OBTAIN  THE  E I GANVALUES  (El) 

C  AT  EACH  XI-LEVEL 

C 

IX=1 

AC=R23 

BC=R43 

CALL  L0SIM(F2I ,H1 I , VI ,TI ) 

eh  d=hii 

I  X  =  2 

S IGMA  C  2  >  =  ( GX 1 13 ( 2) +GX 1 13( 1))*XI(2)/2.0D0 


i  «•  •  i  r  "  -  v‘  I 


\  T i  Q!  !T  J.  J  3 A 
3f  i  JU3JAD  I  .  A  _ 
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idi.FsrSt*  •  >  ta 'prr  £0i 

fl/1 1  .-SI 

mi  ix ,m  >ixd  ircj«?)  3A3.«» 

( r  -i  ‘  ?  i.  -m  i  >  r  x-  (  m  )  i  o 

\  -  ■ 

\  r  ;>  i 

(  4  )  I  X  =  A  ?  AO 

, 

u  '  ;  i  >  o 

£  )ixo=m  f: 

f  30 

(in  Txg*(i:  »ox’ 

.  E  \  t  (  (  -r  )  :  X.  -  (  !  )  I  )  ♦  <  !  -S I  )  J  V0-  1 1  -  U  U  r  0 
».'><•  V  '  -  '  •  •  -  <  .  -  M 

£  -fST  )  IX0*CXI  )£  IX 

-  *  (  i-m  j  ix0=(  i-r  n  utx; 

r  •-•  (  -s '  )  )  •:  '  £  .■ 

ies»i=i  a  on 
.  -  ( '  Mi'  r 

O'  .  =  <  !  M  IMf 

OC  .  0  M  1  )  I  IM1H.5 
<  00.0«U  )I  IMOH2 

no.;  =?a 

iyaj  /a  on  inn  onivire  t  <at2 

a  ,  i>  v  a  !  v  t  -nr.  t-.!  .  v  ”•  i  j*  -in  ■  T-uosf'ua 


f  SflOA 
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m  mn 
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AC=(2.0*(GXI(3)-GXI(  1)  )/GXI  (  2  )  **R43/9.  0D0+R43 
1/GXI13I2) ) *SIGMA(2 ) /X 1(2) / PR+R23/ PR 
CC=R43*SIGMA<  2) /GX 1 13 ( 2 ) /XI { 2 ) /PR 
BC=CC+1 .000/PR 
DC= 1 • CDG-X  1(2) 

EC=AC*2 .ODO 

HC=CC*PR 

GC=BC*PR 

FC=R23*(GXI { 3 )  — GX I ( 1 )  ) *  SI GM A ( 2 ) /GX I (2)**R43 
1/XI (2 ) +  HC 

CALL  LOS  I M ( F2 I»HI I  *  V  I  I  *  T I  I) 

E  I  (  2  )  =  H 1 1 

DO  5  I M=3  *  NOODt  2 

S  S=  SS+ ( GX I  1 3 (  I M— 2 ) +4 • ODO*GX  1 1 3 ( IM-1)+GXU3( IM) ) 
1/3. ODO 

SIGMA  ( IM)=SS*DIF( IM) 

BC=5S^R23/GXI 13 (IM) 

CC= 1 • ODO-X  I  ( IM) 

EC  =  BC=MGXI(  IM  +  1  )— GX  I  (  IM-1 ) ) /GX I ( IM ) 
AC=(EC+9.0D0*8C+2.0D0)/(3.ODO*PR) 

DC= AC*2 .ODO 
I  X=  I  M 

CALL  LOSIMCF2I tHlItVI II ,TII I) 

El ( IM )=H1 I 
I Y=I M+l 

IF  (IY  .EQ.  NM)  GO  TO  5 

ST=SS+(GX  1 1 3  (  IM  )+GX  1 1 3 C  IM  +  1  ))/2.0D0 

SIGMA  II M+l )=ST*OIF{ I M  +  l ) 

BC=ST*R23/GX  1 13 ( IM  +  1 ) 

EC=BC*(GXI  <  IM+2  J-GXI  (  IMU/GXIt  IM  +  1) 
AC=(EC+9.OD0*BC+2.0DO)/(3.0DO*PR) 

CC=l.CDO-XI  ( IM+1 ) 

DC- AC*2.0DO 
I  X=  I M  + 1 

CALL  LCSIM  (F2I»HlI*VIIIf TI II > 

El I  IM+1 )=H1 1 

5  CONTINUE 

IF  (IV  .EQ.  NM)  GO  TO  6 

SL=SS  +  ( GX 1 1 3  C NM-2 ) +4. ODO*GX 1 13 ( NM- 1 ) +GX 1 13( NM) ) 
1/3. ODO 
GC  TO  7 

6  SL=  SS+ ( GX I  1 3 ( NM ) +G X I 13(NM-1 )) /2.0D0 

7  SIGMA(NM)=SL*OIF(NM) 

BC=  SL*R23/GX I  13 (NM) 

EC=BC*(3.0D0~(4.0D0*GXI (NM-l)-GXI (NM-2 ) )/GXI(NM) ) 
AC=( EC+9.0D0*BC+2.0D0 )/( 3.0D0*PR) 

CC=1. ODO-X  I  (NM) 

OC=  2 • ODO* AC 
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I  X  =  N  M 

CALL  LCSIM  (F2  I  ,H1  I , VI I  I, TI 1 1) 

EIINM )  =  H1  I 

SOLVE  THE  CONDUCT  ION  PROBLEM  WITH  A  FOURTH-ORDER 
RUNGE-KUTTA  ROUTINE  AND  A  SHOOTING  METHOD.  START 
THE  INTEGRATION  FROM  THE  ROOT  OF  THE  FIN  WITH  THE 
BOUNDARY  CONDITION  GXI(NM)=1.0  AND  A  GUESSED  VALUE 
OF  THE  TOTAL  HEAT  TRANSFER , HE  AT f NM ), OF  THE 
OVERALL  PROBLEM 

GDHNM)  =  1  .ODO 

8  HEAT ( NM )=2 .ODO 

DIF  SUV  I  1 )  =  (4. ODO*SIGMA( 2) **0.75D0-S IGMAI3) 

1**0  .7 5D0 ) /2.OD0/XI ( 2) 

DIFS  UM{2)=SIGMA{3)**0. 7500/2*0  DO /XH2  ) 

DO  9  1=3, MN 

DIFSUMt I)=I SIGMA U +  1 )**0. 75D0-S IGMA { I-L)**0. 75DO) 
1/2*  ODO/X 1(2) 

9  CONTINUE 

DIFSUM{NM)={3.0D0*SIGMACNM) **0. 75D0-4.CD0 
l*SIGMA(NM-l)**0.75D0+SIGMAf  NM-2 ) **0 . 75 DO ) 
2/2.0DC/XU2) 

10  00  11  I  X  =  2  , NM 
I  R=NM  +  2- 1  X 

GDI  1(  IR  )  =  HE  AT  {  IRJ/XK  IR) 

CL=CN*E I ( IR)*DIFSUM( IR J*GAMA43**0 .75D0 
GOO  1 1=RAT ( IR)*(G0I1< IR ) +CL*GD  It IR) ) 

GDK11  =  D IF (  IR ) *CL*GDI { IR) 

GDK 12=-D I F { IR)*GDI1{ IR) 

GDHF=GD I { IRI+GDKI2/3.000 
GDHH1  =  GD 1 1 (  IR  )+GD01 1/3. ODO 

CL=CN*GXK13 ( IR-1 ) * ( E I (I R ) -{ E I U R ) -E I (  I R- 1 ) ) 

1/3. ODO )/(GAMA3*{ SIGMA! I R- 1 ) +( GX I  13 (  IR-1) +4. ODO 
2*GX  J  13 (  IR-1  J+GXK13  (  IR-1 )  )  *DIFUR)/9.0D0)  )  **R14 
GD021=DIF  {  IR)  *( GCHH1+CL*GDHH  )/ ( XI  UR)-DIF(IR) 

1/3. ODO) 

GDK2 1 =DIF { IR)*CL*GDHH 
GDK  22=— DI F ( I R ) *GDHH1 
GDHh=GD I (  IRH-GDK22-GDK12/3.  ODO 
GDHH1 =GCI 1 (  IR )  +  GD02 1  — GD01 1/3. ODO 
CL=CN*GXJ 131 I  R-  1 )  *  (  E I { I R- 1) +{ E I (IR )-E I ( IR-1 ) ) 
1/3.0D0)/(GAMA3*(SIGMA< I R- 1 ) +{ GX J 13 ( IR-1) 

2+GX  1 1 3 <  IR-1 )  )  *DIFUR)/6.0D0  )  ) **R14 
GD03 1  =  DIF {  IR)*(GCHHl+CL*GDHH)/( XI  UR) -DIF  I IR)*R23) 
GDK31=DIF< IR)*CL*GDHH 
GDK  32=— D I F ( IR )*GDHH1 
GDHH=GD  I  ( IR ) -GDK 2 2  +GDK32+GDK1 2 
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GDHH1=GCI 1 (  IR  J-GD021 +GD03 i+GDOl 1 
CL  =  CN  *E I (  IR-i)*OIFSUM(  I R-  1 )  *G  AM  A43**0 . 75D0 
G  CK4 1  =  0  I F ( IR ) *CL*GDHH 
GDK42=-CIF( IR)*GDHH1 

GDI (  IR- 1 )  =  GD I <IR)+(GDK12+3.0D0*(GDK22+GDK32 ) 
1+GDK42 ) /8  *000 

HE ATI  IR-1  ) =HE AT l IR )  +  { G0K1 1  +  3 .0D0* ( GDK2 I+G0K3 1 ) 
1+GDK41 )/8.0D0 

11  CONTINUE 

CHECK  THE  SATISFACTION  OF  THE  BOUNDARY  CONDITION 
AT  FIN  TIP, HE  AT  II )=Q . 0, AND  THE  CONVERGENCE  OF  GXI 

D I  =  HE  AT { 1 ) 

IF  (01  .LT.  0.0  DO)  GO  TO  19 
IF  (DABS(DI)  .GT.  1. 00-51  GO  TO  19 
DO  12  1=1, MN 
RES  I=GXI (  D-GDUI) 

IF  (DABS(RESI)  .LT.  l.OD-5)  GO  TO  12 
I  J  =  2 

12  CONTINUE 

IF  ( I J  .GT.  1)  GO  TO  24 

13  WRITE  (6, 104) 

104  FORMAT! 10 X, * GXI (IX) • ,14X,*XI{ IX)',16X, '  E I  (  I  X  ) ', 
113X,  *  IX*  ,  I2X,  *  SIGMA!  IX)  •  ,12X,»  GXOUX)*) 

DC  14  I X= 1 , NM 

WRITE  (6,105)  GXI ( I  X) , XII IX) , El ( I  X) , IX  ,SIGMA(IX)  , 
1GX0<  IX) 

105  FORMAT  (■  *,3D23.16,  14,  2D23.I6) 

14  CONTINUE 

WRITE  (6,106)  ITERA,LCOP 

106  FORMAT  ('  I TERATION=  *,I6,*  LOOP=  *,I4) 

CHECK  THE  MATCHING  OF  INTERFACIAL  TEMPERATURE 

DC  15  1=1, MN 
RES  I=GX I (  I)-GXO(I) 

IF  (OABS(RESI)  .LT.  l.OD-5)  GO  TO  15 
I  K=2 

15  CONTINUE 

DO  18  I X=1 » NM 
GX  =  GX  I  ( IX ) 

IF  (GX  .GT.  l.OD-18)  GO  TO  17 

GXI ( I  X )  =  1 • OD- 18 

WRITE  (7,103)  GXI ( IX),  XI(IX) 

GO  TO  18 
17  GXI ( I  X )  =  GX 

WRITE  (7,103)  GX I (  IX),  XI(IX) 
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18  CONTINUE 
I  K=  1 

LQOP=LOCP+ 1 


IF 

(IK  .GT. 

1  ) 

GO 

TO 

KC  = 

KO  +  1 

IF 

(KO  .EQ. 

2) 

GO 

TO 

WRITE  (6,107) 

107  FORMAT!'  PROBLEM  IS  SOLVED  ') 

STOP 

19  I TERA= I TERA+1 

OF  THE  TOTAL  HEAT  TRANSFER, HE AF (NM ), SO  THAT  THE 
THE  TOTAL  HEAT  TRANSFER  R AT E, HE AF ( NM ), SO  THAT  THE 
SHOOTING  METHOD  FOR  MODIFYING  THIS  El GANVALUE  CAN 
8E  APPLIED 

GDI  I ( NM  )=  1 *000 
HEAF(NM)=HEAT(NM)*1,1D0 

20  DO  21  I  X=2, NM 
I R  =  N  M+2- 1 X 

GDI  Il{  I R  )  =HEAF<  IRJ/XIUR) 

CL=CN*E  I  {  IR  )*DIFSUM(IR  )*GAMA43**0. 75D0 
GD011=RAT ( I R ) * ( GD I  Il(  IR  ) *CL  *GD I  I (  IR ) ) 

GDKU=DIF(  IR  )  *CL*GD  1 1  UR) 

GDK 12=-CI F(IR)*GDII1( IR) 

G  DH  H=GD I  I ( I R ) +GDK1 2/3 .000 

G0HH1=G0II 1( IR) +GD0 11/3. ODO 

CL=CN*GXK13( IR-1 )*(€! ( IR )-( El ( IR)-E  I(  IR-1) ) 

1/3. 0D0)/(GAMA3* ( SIGMA ( IR-1 M-( GX 113 ( IR-1 ) +4. ODO 
2*GXJ13( IR— 1 J+GXK13I IR-1) )*DIF( IRI/9.0D0) )**R14 
GD02 1 =D I F {  IR) *( GDHHl+CL*GDHH)/( XI  (IR)-DIF(IR) 

1/3. ODC) 

GDK21=DIF( IR ) *CL*GDHH 

GDK22=-CI F ( I R  >  *GDHHI 

GDHH=GD 1 1 (  IR) +GDK22-GDK 12/3 .ODO 

GDHH1=GD I I 1( IR ) +GD02 1-GDO 1 1 /3 . ODO 

CL=CN*GXJ13( I R-l )* (FI  (  I R- 1 ) + ( E I( I R )-E K IR- 1 ) ) 

1/3. ODO )/(GAMA3*( SIGMA! IR-1) +GXJ13( IR-1 ) 

2+GX I 1 3 ( IR-1))*DIF( IR)/6.0D0) )**R14 
GD031  =  DIF  (  I R )  *  (  GDHH1  +  CL #GDHH  )/ ( X  I (  IR)-CIF(  IR)*R23) 
GDK3 1=DIF ( IR ) ^CL^GDHH 
GDK32=— Cl F ( IR ) *GDHH1 
GDHH=GDI I ( IR ) -GDK22  +GDK32  +GDK12 
GDHH1  =  GD  III  (  IR)-GD021+GD0  3U-GD011 
CL=CN*EI ( IR-1 )*DIFSUM( IR— 1) *GAMA43**0. 7500 
GDK41 =D IF ( IR ) *CL*GDHH 
GDK42=-DI F ( IR ) *GDHH1 

GCI  II  IR-1 )  =  GD  I  I ( IR  )+( GDK 12+3. ODO* ( GDK2  2+GOK32 ) 
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1+GDK42) /8.0D0 

HEAF(  IR-1 }=HEAF( IR)-MGDK1 1+3.000* { G0K2  1+GDK31 ) 

1 +GDK4 1 ) /8 . 0D0 

21  CONTINUE 

CHECK  THE  SATISFACTION  OF  THE  BOUNDARY  CONDITION 
AT  FIN  TIP ,HE AF m =0. 0, AND  THE  CONVERGENCE  OF  GXI 

D I  I =HEAF ( 1  ) 

IF  IDII  .LT.  0,000)  GO  TO  37 
IF  (DABS(DII)  .GT.  1.0D-05)  GO  TO  37 

22  DC  23  I X= 1 » MN 

RES  1=GX  I  (  I  X )-GD I  I {  IX) 

IF  (DABS(RESI)  .LT.  1.00-05)  GO  TO  23 
1 1=2 

23  CONTINUE 

IF  (  IL  .GT.  1 )  GO  TO  29 

FEEC  THE  CONVERGED  GXI  OF  CONDUCTION  EQUATION  BACK 
TO  THE  BOUNDARY  LAYER  PROBLEM 

GO  TO  13 

IF  GXI  DOES  NOT  CON VERGE , GE NERATE  A  NEW  GXI  AND 
CALCULATE  THE  COEFFICIENTS  FOR  THE  NEW  ITERATION 
OF  CONDUCT  I CN  LOOP 

24  DO  27  IX=1,MN 

GX=  <GD I ( I  X ) *GX I (IX) ) /2 . ODD 
IF  (GX  .GE.  l.OD-18)  GO  TO  25 
GXI { IX)=1.0D-18 
GO  TO  26 

25  GX I  I  I  X )=GX 

26  GX 1 13  < IX ) =GX I ( I  X ) **R13 

27  CONTINUE 

DO  28  IX=2,NM 

GXJ ( IX— 1 )=GXI ( iX-l)+(CXI( IXJ-GXII IX-1 ) )/3.0D0 
GXK I IX-1 )=GXJ ( I X-l )  +  <  GXI ( IX  I— GXI < IX-1) )/3.0D0 
GXJ  131  IX-1  )=GXJ  (  IX-1)^R13 
GXK13 ( IX-1 )=GXK (IX-1) ^*R13 

28  CONTINUE 
I  J=  1 

GO  TO  34 

29  DO  32  I  X= 1 »  MN 

GX= ( GO  1 1 ( I  X ) +GX I C IX) )/2.0D0 
IF  (GX  .GE.  l.OD-18)  GO  TO  30 
GXI ( I  X )  =  1 . OD-  1 8 
GC  TO  31 
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30  GXI (I  X  )=GX 

31  GXU3I  IX}=GXI(  I  X)**R13 

32  CONTINUE 

DO  33  IX=2,NM 

GXJ  {  IX-U  =  GXI  II  X-l  )+IGXII  IXI-GXI  {  IX-L  )  )/3.0D0 
GXKUX-1  )=GXJ{  I  X-l  )  +  I GX  I  (  IXI-GXII  IX-1 >)/3. ODO 
GXJ  13  {  IX-1  )=GXJ  (IX-1)  **R13 
GXK13  C IX-1 )=GXK< IX-1) **R13 

33  CONTINUE 
I  L=  1 

34  S IGMA 1 2 )= I GX 1 13 { 2 ) +GX I  1 3 C  1)  ) *XI { 2 ) /2. ODO 
DO  35  1=3  *  NOOD»  2 

S IGMA ( I J  =  S I GM  A I  I-2J  +  IGXI 13(  1-2 ) +4. ODO+GXI 13 ( I - 1 ) 

1  +  GX 113 ( U ) / 3 • 0D0*D  IF  {  I  ) 

SIGMA (I+l)=SIGMA(I)+{ GX  1 1 3 (  l )+GXI13t I+1))/2.0DO 
1*D I  F  (  1  +  1) 

35  CONTINUE 

IF  II V  .EQ.  NM)  GO  TO  36 

SIGMA  (NM)=SIGMA(NM-2)  +  1  GX  II  3(  NM-2  1  +  4.  ODO 
1*GX  1131  NM-1H-GX  113  I  NM  ) ) /3  .0D0*DIF < NM  ) 

36  I TERA=I TERA+1 
GO  TO  8 

37  DI  =  HEAF (1 ) -HE  AT ( 1 ) 

IF  IDABSIDI)  .LT.  1.0D-20)  GO  TO  22 

CALCULATE  AN  IMPROVED  E I G ANVALUE , HE AF ( NM ) , B Y  MEANS 
OF  THE  SHOOTING  METHOD 

GDCGR-HEAT ( NM )  +  I HEAT { NM )-HEAF (NM) )/(HEAF( 1) 

1— HEAT I i)) *HEAT{ 1) 

HEAT ( NM )= HEAP (NM ) 

H  EAF ( NM )  =  GDCOR 
HE  ATI  1 ) =HE  AF ( 1 ) 

I TERA= I  TER A  +  l 
GO  TO  20 
END 
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C 

c 

C  SUBROUTINE  PROGRAM  FOR  INTEGRATING  THE  BOUNDARY 

C  LAYER  EQUATIONS  WITH  "LOCAL  SIMILARITY”  TECHNIQUE, 

C  THE  INTEGRATION  IS  ACCOMPLISHED  USING  A  FOURTH- 

C  ORDER  RUNGE-KUTTA  ROUTINE  AND  THE  CONVERGENCE 

C  IS  BASED  ON  THE  USE  OF  A  LEAST-SQUARES  METHOD 

C 
C 

SUBROUTINE  LOSIM( F2 I , HI  I , FUNV , FUNT ) 

REALMS  SFOMI 1(251) , SF1MI II 2 51) , SH0MI 1  I  251 ) , 

ISHIMI II 251 ), SFOMI 21 251) ,SF1M12 I  251) ,SHOMI2I  251)  , 
2SH1MI2I251)  ,F2I251) , FI (251) , F I  251 ) , HI  I  251 ) , 

3H(25i ) ,SF2< 251) ,SF I( 251) , SF{251) , SHI I  2  51) , SHI  251) 
REAL*8  F2MI4) , F  1M I  4 ) , HIM I  4) ,  HM  (  4)  ,  A  21  4  )  ,  A 1  {  4)  , 

1GXI I  ICO) ,XII100),AC,BC,CC,DC,EC,FC,GC,HC,PR, 

2ERRDI F, AERRDI 

REAL* 8  F2I ,Fl I, FI, HI I ,HI , FF2,FF1, FF,HH1,HH,T, V, 
1K11,K12,K13,K14,K15 ,K21,K22  ,K23,K24  ,K25,K31 ,K32, 
2K33 ,K34,K35,K41 , K4 2 , K43 ,K44 ,K45 , ERROR, ERRORO, 

3 ERR, ERR AT , DEN  ,  DF2 , DH1 , ET A ,G , T EMP, VELO 
REAL  C,P,Q»R»U»  W,YS 

COMMON  SF2,SF1, SF , SHI ,SH, SFOMI 1,SF1MI 1 ,SH0MI1 , 
ISHIMI 1, SFOMI 2 ,SF1MI2,SH0MI2,SH1MI2, PR,  AC, BC,CC, 
2DC, EC  »  FC , GC , HC  »  GXI , XI , I  X, I , J ,KO , KOA , KGB , KOC , 
3KCD,K0E 
G  =  0*  0400 
ERRORO=O.ODO 
M-5  1 

H 1= 1, ODO 
F 1=0, CDO 
F 1 1=0 ,000 

50  J  =  1 
N=M— 1 

START  INTEGRATING  THE  BOUNDARY  LAYER  EQUATIONS 

F2I  1 ) =F2 I 
F  1  (  1 )  =F 1 1 
FI  1)  =  FI 
H1I  1)=H1I 
HI  1)=HI 
GO  TO  53 

51  IF  IJ  . GT ,  2)  GO  TO  52 

START  INTEGRATING  THE  PERTURBATION  EQUATIONS  OF 
THE  BOUNDARY  LAYER  EQUATIONS  {DUE  TO  THE  VARIATION 
OF  EIGANVALUE  F2I) 
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F2C1I=1.0D0 
Fl( 13  =0.000 
F ( 1)=O.ODO 
HU  1  )=0.0D0 
H  {  1  )=0. 0D0 
GO  to  53 

52  CONTINUE 

START  INTEGRATING  THE  PERTURBATION  EQUATIONS  OF 
THE  BOUNDARY  LAYER  EQUATIONS  (DUE  TO  THE  VARIATION 
OF  E I GANVALUE  H1I) 

F2( 11=0.000 
F 1 ( 1 ) =0. 0D0 
F { 1 )=  0. 0D0 
HI ( 1 ) =1 .0D0 
H ( 1 ) =G*  0D0 

INTEGRATE  THE  EQUATION  NITH  A  FOURTH-ORDER 
RUNGE-KUTTA  METHOD 

53  DO  54  1=1, N 
FF2=F2< I > 

FF1=F 1( I ) 

F F=F (  I) 

HH=H (  I  ) 

HH1=H1( I ) 

CALL  FUNV(V,FF2»FF1»FF,HH1,HH) 

CALL  FUNT (T ,FF2,FF1 ,FF,HH1,HH) 

K 1 1=G^V 
K 12=G*FF2 
K13=G*FF1 
K14=G*T 
K15=G*HH1 

FF2  =  F  2  (  I  H-Kll/3  .000 

FF1=F 1 ( I ) +K12/3 .ODO 

F  F  =  F (  II +K 13/3. ODO 

HH=H { I ) +K15/3 .000 

HH1=H 1 ( I ) +K14/3 .000 

CALL  FUNV ( V  ,FF2 ,FF 1 ,FF ,HH 1, HH) 

CALL  FUNT  (  T  ,  F F2. »  F F 1 ,  F F ,  HH1 ,  HH  ) 

K21=G*V 
K  22  =  G*FF2 
K23=G*FF1 
K  24=G*T 
K25=G*HH1 

FF2=F2 ( I J+K21-K 11/3 .000 
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FF1=F  H I ) +K22-K 12/3 .000 
FF=F (  I  3+K23— K 13 /3 *  0D0 
HH=H{ I ) +K25-K15/3 .  ODO 
HH 1=H 1 ( I ) +K24-K14/3.0D0 
CALL  FUNVI V,FF2,FF1,FF,HH1,HH) 

CALL  FUNT ( T, FF2 ,FF1 , FF, HH1,  HH) 

K31=G*V 
K32=G*FF2 
K33=G*FF1 
K  34=G*T 
K35=G*HH1 

FF2=F2U)-K21+K31+K11 

FF 1  =  F  1  (  I  )— K22+K32+K 12 

F  F  =  F (  I ) -K2  3+K33  +K1 3 

HH=H ( I J-K2  5+K35+K1 5 

HH 1=H 1 ( I ) -K24+K34+K 14 

CALL  FUNV (V»FF2*FFlfFF»HHlfHH) 

CALL  FUNT(T,FF2,FFl,FFfHHl,HH) 

K41=G*V 
K42=G*FF2 
K43=G*FF 1 
K44=G*T 
K45=G*HH1 

F2( I + 1 ) =F2 I I) +IK11+3.0D0* IK21+K31 )+K41 )/8.0D0 
Fit  I  +  1)  =  F1(  I )  +  (K12  +  3.0D0’MK22+K32  )  +  K42  )  /  8 . ODO 
FU+n=FUH(  K1 3+3 . ODO* I K23  +K33  )+K4 3)78.0 DO 
HHI+  1)=HH  I)  +IK14+3.0D0*  ( K24+K 34 )+K44 ) 78.000 
H  1 1  +1  )  =  H  (  I  )  +(  Ki  5+3  •  ODO*  (K25+K35J+K45)  /8.0D0 

54  CONTINUE 

IF  (J  .GT.  1)  GO  TO  56 
DO  55  1=1, M 
SF2 ( I )=F2 ( I ) 

SF1 ( I  )  =  F1 { I ) 

SFC  I)=F(I) 

SHI  I ) =H I  I ) 

SHI  I  I  )  =  H1 I  I ) 

55  CONTINUE 

56  F1MIJ )=F1 (M) 

F2M  (  J )=F2 I M ) 

H  f*  (  J)  =H(M) 

H1MI J  )  =  H1 1 M) 

J  =  J  +  1 

IF  (J  .GT.  3)  GO  TO  57 
GO  TO  51 

CHECK  THE  MEAN  SQUARE  ERROR 


57  0=F  1M I  1 )**2 
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1  +  11 
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P  =  F2M 1 1 }#*2 
Q  =  HM(  l  )**2 
R=H1M( 1 )**2 
ERRGR=G+P+Q+R 

IF  (ERROR  .LT.  1.0D-I0)  GO  TO  58 
ERR=( ERROR-ERRORC) /ERROR 
ERRAT=DABS (ERR) 

GO  TO  59 

58  ERR AT=0  *000 
FRR=0. 000 

59  ERRDI F=ERROR-ERRORO 
AERRCI=DABS( ERROIF ) 


I 

F 

( M 

•  EQ. 

101) 

GO 

TO 

61 

I 

F 

(M 

•  EQ. 

151) 

GO 

TO 

63 

I 

F 

(M 

.FQ. 

201) 

GO 

TO 

65 

I 

F 

(M 

.EQ. 

251) 

GO 

TO 

67 

I 

F 

(ERROR  , 

.LT. 

1.00 

-05)  GO 

TO  60 

I 

F 

(AERROI 

.LT. 

1.00- 

05)  GO 

TO  60 

IF  (ERRAT  •  GT  .  1.00-05)  GG  TO  68 

60  ERRGR0=0. 000 
M  =  1 01 

GO  TO  69 

61  IF  (ERROR  .LT.  1. 00-05)  GO  TO  62 
IF  (AERRDI  .LT.  1.0D-05)  GO  TO  62 
IF  (ERRAT  .GT.  1.00-05)  GO  TO  68 

62  ERR OR 0=0. 000 
M=  15 1 

GO  TO  69 

63  IF  (ERROR  .LT  .  1.00-05)  GO  TO  64 
IF  (AERROI  .LT.  1. 00-05)  GO  TO  64 
IF  (ERRAT  .GT.  1.00-05)  GO  TO  68 

64  ERROR 0=0. 0 DO 
M=201 

GO  TO  69 

65  IF  (ERROR  .LT.  1.00-05)  GO  TO  66 
IF  (AERRDI  .LT.  1.00-05)  GO  TO  66 
IF  (ERRAT  .GT.  1.00-05)  GO  TO  68 

66  ERR0RG=0. 0 DO 
M  =  251 

GO  TO  69 

67  IF  (ERRAT  .LT.  1.0D-10)  GO  TO  70 
IF  (ERROR  .LT.  1.00-10)  GO  TO  70 
IF  (AERRDI  .LT.  1.00-10)  GO  TO  70 

68  ERRGRC=ERROR 

GENERATE  AN  IMPROVED  SET  OF  EIGANVALUES  F2I  AND 
H  1 1 
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69  Alt  1)  =  F1M(  1)*F1M(  2)  +HM{  1 )  *HM(  2)  +F  2M  t  1)  *F2M(  2  ) 
l+HIMt 1)*H1M<2) 

A1(2)=F1M(2)**2+HM(2) **2+F2M( 2) **2+Hl M ( 2 ) 

A1  (  3  )  =  F 1 M  (  2  )  *  F 1  M  (  3  )  +  HM  (  2  )  *HM(  3)+F2M(  2) *F2M( 3) 

1+H 1 M ( 2 ) *H 1  M  ( 3  ) 

A 2  t  1 )  =  F 1 M  (  1)*F1M(  3M-HM(  1)  *HM(  3)  *F2M(  1)*F2M(3) 
l  +  HIMt  1)*H1M<3 ) 

A2  t  2 ) = A 1 ( 3 ) 

A2  (  3 )=F 1M (  3)**2+HM(  3 )  **2+F2M(  3  )  **2  +  H  1 M  1 3  )  **  2 
DEN  =  AU2)*A2(3)-A1  (3)  *A2(2) 

OF  2=A 1 1 1)*A2(3)-A1(3)*A2( 1) 

DH 1=A 1 1  2 ) * A2 1 1 ) —A 1(11 *A2 ( 2 J 
F2  I=F2 I -OF2/DEN 
H  l  I=H 1 1-DH 1 /OEN 
GG  TO  50 

70  F  2  I =SF2  ( 1 ) 

H 1  I  =  SH1 ( 1) 

IF  THE  CONVERGED  SOLUTION  OF  THE  OVERALL  PROBLEM 
IS  OBTAINED, PUNCH  THE  DATA  OF  THE  BOUNDARY  LAYER 
PROFILES  AT  THOSE  SPECIFIED  XI-LEVELS  CF  KOA, KOB, 
KOC ,KQD , AND  KOE 

DC  72  1=1, M, 2 
I  E=  I  —  1 

ET  A=G*CFLQAT  t IE ) 

IF  (IX  .EQ.  KOA)  GO  TO  71 

IF  (IX  .EQ.  KOB )  GO  TO  71 

IF  (IX  .EQ.  KOC )  GO  TO  71 

IF  (IX  .EQ.  KOD)  GO  TO  71 

IF  (IX  .EQ.  KOE)  GO  TO  71 

GO  TO  72 

71  IF  (KO  .EQ.  1)  GO  TO  73 

WRITE  (6,200)  IX,GXI< IX) , XI ( IX) 

200  FORMAT  ('  IX=',I4,*  GX  1=  * ,  D23.  3.6 ,  1  XI=*  D23.16) 

WRITE  (6,203) 

203  FORMAT ( * 0  ET A* , 8X,  *F  • , 14X, *F 1 • , 13X  , 

1 ' F 2 1  , 13X , *  H  * , 14X, *  HI  *  ) 

T  EMP=  SH ( I  ) 

VEL0=SF1 ( I ) 

W=SNGL( TEMP) 

U=SNGL ( VFLO ) 

Y  S=SNGL ( ETA ) 

WRITE  (7,204)  L,W,YS 

204  FORMAT  ( '  *  ,3E13.6) 

WRITE  (6,205)  ETA  ,  SF ( I) , SF 1 ( I ) , SF2 ( I ) , SH ( I ) , 

1SH1 ( I ) 

205  FORMAT  (»  ' , 4X , F10 . 4, 5E 15 .3 ) 
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72  CONTINUE 

WRITE  <6,206)  ERROR, ERRAT ,AERRDI 
206  FORMAT <  *  ERRGR= • , 023.  16 , *  ERRAT=  *,023.16, 
1*  AERRDI=  *,023.16) 

73  DO  74  1=1, M 
SF1MI2C T )=SF1 MI 1 < I ) 

SF0MI2U  )  =  SFOM 1 1  <  I  ) 

SH1MI 2< I )=SH1MI 1< I ) 

SH0MI2 ( I )=SHOMI 1 ( I ) 

SF1MIU  I )  =  SF1  ( I  ) 

SF0MI1U  )  =  SF  (  I) 

SH1MIK  I  )=SH1  (  I  ) 

SHOMIli I )  =  SH (  I) 

74  CONTINUE 
RETURN 
ENO 
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FUNCTION  SUBROUTINE  OF  EQUATION  OF  MOTION  AND  ITS 
PERTURBATION  EQUATI ONC J>I )  AT  THE  FIRST  XI-LEVEL 


SUBROUTINE  VIC V , FF2 , F FI , FF , HH 1 , HH ) 

REALMS  FF2,FF1,FF,HH1 ,HH, V 

REAL* 8  SF2(251)  ,SF 1(251  ),  SFC25I  )9 SHU  251)  ,SHC25I)  , 
1SFOMI  1(  251)  tSFIMI  1(  25  1)  ,  SHOMI  1(  251)  ,  SH  IMI  1(  251 )  , 
2SF0MI2I  251)  ,SF1MI2(  25  1)  ,  SHOMI  2 (  25  1)  ,  SH  IMI  2 1  25  L  )  , 
3GXK ICO) ,XI (100),AC,BC,CC,DC,EC,FC»GC,HC,PR 
COMMON  SF2,SF1, SF,SH1,SH, SFOMI 1, SF1 MI  1 , SHOMI 1 , 
1SH1MI1,SF0MI2 ,SF1MI2,SH0MI2,SH1MI2,  PR, AC, BC,CC, 
2DC,EC,FC,GC,HC,GXI ,XI , I  X , I , J  ,KO , KOA , KCB , KOC , 
3KCD,KOE 

IF  (J  .GT.  1)  GO  TO  81 

V^( AC*FF1*FF1-FF2*FF) /PR-HH 

RETURN 

81  V={ BC*FF1*SF1 ( I }-FF2*SF( I )-FF*SF2 ( I ) }/ PR-HH 
RETURN 
END 
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FUNCTION  SUBROUTINE  OF  ENERGY  EQUATION  AND  ITS 
PERTURBATION  EQUAT ION C J>I }  AT  THE  FIRST  XI-LEVEL 


SUBROUTINE  T I ( T , FF2 , F FI ,F F ,HH 1, HH ) 

REAL*8  FF2 ,FF1, FF , HHI ,HH, T 

REALMS  SF2(  251)  VSF1(25II«SF  (251)f  SHK251)  «SH(251)  t 
1SFGMI 1(251) ,SF1 MI  1(251) ,SH0MI1( 251)  ,SH1MI 1( 251)  , 
2SFOMI 2(251) ,SF1M 12 (251) , SHOHI2 ( 25 1 ) ,SH 1MI2(  251 )  , 
3GX I ( ICO) ,XI { 100) , AC,BC,CC  ,DC,EC ,FC,GC ,HC,PR 
COMMON  SF2,SF1, SF , SHI , SH, SFOMI 1 , SF 1  MI 1,SH0MI1, 
1SH1MI 1, SFOMI2 ,SF1MI2, SHOMI2 , SHI  MI  2 , PR , AC, BC,CC, 
2DC, EC  »FC, GCtHCt  GX I , XI , I X , I , J  ,KO , KOA , KCB, KOC , 

3KGD , KOE 

IF  { J  .GT.  1)  GO  TO  82 

T  =-FF^HHl 

RETURN 

82  T=-FF*SH1( I)-HH1*SF(I) 

RETURN 

END 
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FUNCTION  SUBROUTINE  OF  EQUATION  OF  MOTION  AND  ITS 
PERTURBATION  EQUAT ION( J>I )  AT  THE  SECOND  XI-LEVEL 


SUBROUTINE  V 1 1 i V , FF2 , FF1 ,F F , HH1 , HH ) 

REALMS  FF2.FF1, FFfHHl,HH, V 

REAL* 8  SF2I251)  , SF 1 ( 2 51 ) , SF ( 251 ) , SH 1 { 25 1 ) , S H{ 25 1 ) , 
1SFGMI 1C  251) f SF1MI 1( 251) ,SHQMI 1(  251)  ,SH1MI 1( 251) , 
2SFCMI2C251)  ,SF1MI2(  251)  , SHOW  12 (  251)  ,SH1MI2(  251)  , 
2GXI (ICO) ,XI( 100) ,AC,BC,CC,OC,EC,FC,GC, HC, PR 
COMMON  SF2,  SFi,  SF  ,  SHI  ,SH,  SF  OMI  1 ,  SF  1  MI  l,SHOMIi, 
1SH1MI 1, SFOMI2 , SFI MI  2, SHOMI2 , SH1MI 2, PR, AC, BC,CC, 
2DC, EC,FC, GC,HC ,GXI , XI , I  X , I , J ,KO , KOA,KOBf KOC , 
3KQD,KOE 

IF  (J  .GT.  1)  GO  TO  83 

V=AC*FF1*FF1-BC*FF*FF2+CC*( SFOMIK  I  )*FF2 
1-SF1MIK  I  )  *  FF 1 )  -DC  *HH 
RETURN 

8  3  V=EC*FF1*SF1(  I  )-BC*(FF* SF2C  I ) +FF2*SF(  I  )  ) +CC 
1*(  SFOMIK  I)*FF2-SF1MI  1<I)*FF1  )-DC*HH 
RETURN 
END 


:  ■  V  -IX  n/TK>  !?  T A  (  l  'l  >  i'T.  ’  7  ' :  '  T< 


{  .h  ,  J  !,  -1  ,  i  ?  ,  ,V)  J  IV 

( r  i  tib  ,  :i  '  *  , 

.  i  C--  )!.,,!!  )  (i  .<  f -‘S'  ,  (  •*■• )  :  <•<?.  ( !<*.  1. ) 

,  ( j  )  f  r  t  (  r  <  !  ]  1*  i  !  1 1  M  ’  :  ■  ■  t  (  )  :  I 

,  ues)  )  •  •  '  i  '  .  1MCH21 

f'  ,3H,  J,  .  >3,  ,  3 c f  ,  ■  ♦  r  (  ’  l )  !  >  i  t  vr  )  I >  > 

♦  I'  H  ,  i  I'  I  V  ,  {  ;  •  ,H2,  ,  1?  ,  I  2  r  ' 

*  f  •  !  .  -  '  t  .  i  •  '  .r  •  : 

♦  /  tao>t/v  ,  i  ,x  i,  ix,  ix  »oh  *  *o « r>- .  •  ,oos 

■:  >  t  00  XE 

r8  i  T  )0  (I  .1  .  U  ;  I 

i  -f 

hi  *oo-(  i  =i*{  nir  v.iqg-1 

10  +  < <  t  )  ?  ■  -n  ■  '  •  ?  )  ?  -i  )H?*  ri  -v  8 

HH  •  -<  •  •  (  !  )i  I  i  i  ~  -S'  i*{  I  ),  1  I 


■  ■  >  n  .  o 


o  o  o  r>  o  o 


102 


FUNCTION  SUBROUTINE  OF  ENERGY  EQUATION  AND  ITS 
PERTURBATION  EQUAT ION ( J>1 )  AT  THE  SECOND  XI-LEVEL 


SUBROUTINE  T I  I  { T , FE2 , FF1 , F F , HH 1, HH ) 

REALMS  FF2tFFI,FF,HHl,HH,T 

REALMS  SF2 ( 251 ) «SF 1 ( 251 ) «  SF (251 ) » SH1<  251) ,SH( 251 ) , 
ISFOMI  1(251)  VSF1HU(251)  »SHOMI  1(251 )  fSHIMI  1(251)  v 
2SFOM I 2( 251) ,SF1MI2( 251) ,SHOMI2( 251) ,SH1MI2( 251) , 
3GX I (ICO) t XI (100),AC,BC,CC,DC,EC,FC,GC,HC,PR 
COMMON  SF2.SF1, SF , SHI ,SH, SFCMI1 ,SF1MI l,SHOMIl, 

1  SHIM I l,SFOHI2 ,SF1MI2, SHOMI 2 , SHI  MI  2 , PR , AC,BC,CC, 
2DC, EC,FC,GC,HC,GXI ,XI ,  IX, I , J , KO , KOA , KOB ,KOC , 
3K0D,KQE 

IF  (J  .GT.  L)  GO  TO  84 

T=FC*FF1*HH-GC*FF*HH1  ♦HC*  (  S  FQM  I  1  (  I  )  *HH1 
1-SHOMIK  I  )*FF  1) 

RETURN 

84  T=FC*(FF1*SH( I)+HH*SF1 ( I) )-GC*( HH1*SF ( I ) +FF 
1*SH1(  I)  )+HC*(  SFOMIK  I  )*HH  1- SHOM  1 1  (  I  )*FFl) 

RETURN 

END 
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FUNCTION  SUBROUTINE  OF  EQUATION  OF  MOTION  AND  ITS 
PERTURBATION  EQUAT ION ( J>1  )  AT  THE  I-TH  XI-LEVEL 
WHERE  I >2 


SUBROUTINE  VI  I  I ( V ,FF2 ,FF1 ,FF,HH1,HH) 

REAL*  8  FF2 ,FF1,FF,HH1 ,HH, V 

REAL*  8  SF2 ( 25 1),SF  1(251  ),  SF(251),  SH1(  251),  SHI  251)  , 
1SFOMI 1( 251) fSFlMI 1 (251) ,SHO MI  1(251)  ,SH1MI1( 251) , 
2SFQMI 2(251) ,SF1MI2( 251) ,SHOMI2(251) , SH 1MI  2 (  25 1 )  , 
3GX I (100), XI (100), AC,BC,CC,DC, EC  ,  FC,  GC,  HC,  PR 
COMMON  SF2,SF1,SF  ,SH1  ,SH, SF CMI 1 , SF1MI 1 ,SH0MI1, 

1  SHI  MI 1,SF0MI2,SF1MI2»  SHOM  T2,SH1MI 2, PR ,AC,BC ,CC, 
2DC,EC,FC,GC,HC, GXI , XI ,IX,  I , J , KO, KGA ,K0  B»KOC , 

3KCD, KOE 

IF  (J  • GT  .  1)  GO  TO  85 

V=AC*FFl*FFl-( BC*( FF1*(4.0D0*SF1MI 1 ( I )-SF 1M 12 ( I ) ) 

1  +  FF2*(3*ODO*FF-4.0DO*SFOMI1(I )  + SFOM I  2 (  I  ) ) ) 

2+ FF*FF2 )/PR-CC*HH 
RETURN 

8  5  V=DC*  F  F 1*  SF 1 (  I)-( 3C*( FF 1* ( 4.0D0*SF 1MI 1 (I) 

1- SF 1M 12 ( I ) )+FF2*(3»0D0*SF( I  )-4. ODO*SFOM  1 1 ( I  ) 

2+ SFOM  12 ( I  ) )+3. ODO*FF*SF2(  I ) )  +  SF ( I ) *FF2  +  FF*SF2 ( I ) ) 
3/ PR— CC*HH 
RETURN 
END 
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C 

c 

C  FUNCTION  SUBROUTINE  OF  ENERGY  EQUATION  AND  ITS 

C  PERTURBATION  EQUAT ION < J>1 )  AT  THE  I-TH  XI-LEVEL 

C  WHERE  I >2 

C 

c 

SUBROUTINE  T 1 1  I { T , FF2 ,F FI ,F F , HH1 , HH ) 

REAL*  8  FF2,FFI, FF,HHI ,HH, T 

REAL* 8  SF2 C  25 1 ) , SF  1  { 2 51 ) , SF < 25 1 ) , SH 1C  25  I ) , SH ( 25 1 ) , 
1SFOMI 1C  25  1 1 .SF1 MI  1(251) ,SH0MI1 (251)  ,SHiMI 1{ 251) , 
2SFCMI2C  251) ,SF1MI2<  251) ,SHOMI2C  251) ,SH1MI2( 251) , 
3GXIC ICO  )  ,XI (100 ) , AC,BC,CC,D€, EC,FC,GC, HC,  PR 
COMMON  SF2, SF1, SF  ,SH1 .SH, SFGMI 1. SF1MI 1 ,SH0MI1, 
1SH1MI 1,SF0MI2,SF1MI2, SHOM 12 , SHI  MI  2 , PR ,AC,BC,CC, 
2DC. EC,FC,GC,HC, GXI, XI  ,IX, I, J, KO, KOA,KGB,KOC , 

3KQD , KCE 

IF  ( J  • GT •  1)  GO  TO  86 

T=£C*FF1*HH-FF*HH1+BC*( FF 1* C 3 .0D0*HH-4.0D0 
1*  SHOM  Il(  I  )  +  SHOM  12  (  I  )  )  -HH1*(  3.0D0*FF~4  .ODO 
2* SFOM  1 1 ( l )  +SFQM I2C I  ) )  ) 

RETURN 

86  T=EC*CFF1*SH(  I)  +HH*SF1 {  I)  )-FF*SHl(I  )  -HH  1*SF  (  I  ) 
1+BC*( FF1* ( 3.0D0*SH{  I )-4.OD0*SHGMI 1C  I ) + SHOM I 2( I )  ) 
2+3. ODO*HH*S Fill  )-HH 1* ( 3 .QDO *SF (  I )-4.0DC*SFOMI 1C  I ) 
3+SFOM  12 C I ) )-3. CD0*FF*  SHI ( I ) ) 

RETURN 

END 
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